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CHAPTER  1 
INTRODUCTION 


The  result  of  this  research  is  a continuation  of  the  development  of  the  image  algebra,  an 
algebraic  structure  for  digital  image  processing  applications  [1],  Specifically,  a study  of  the 
use  of  the  image  algebra  as  a mathematical  model  and  tool  for  the  development  of  linear 
and  nonlinear  recursive  image  processing  algorithms  is  presented. 

The  image  algebra  is  an  algebraic  language  for  specifying  image  processing  applications. 
It  represents  an  effort  to  provide  a unified  mathematical  theory  for  image  processing  algo- 
rithm specification,  optimization,  rapid  prototyping,  and  parallel  processing  [2, 3, 4, 5, 6].  The 
development  of  the  image  algebra  is  the  direct  result  of  the  increasing  demand  in  digital 
image  processing  community  for  new  and  more  powerful  algebraic  systems  that  can  provide 
a theoretic  foundation  for  studying  problems  in  image  processing,  express  commonly  used 
image  processing  algorithms,  and  model  parallel  image  processing  architectures. 

1.1  Development  of  Image  Algebra 

The  search  for  high  speed  parallel  computers,  and  parallel  image  processing  architectures 
in  particular,  can  be  traced  back  to  J.  von  Neumann’s  work  on  cellular  automata  in  the 
1950’s  [7].  The  model,  based  on  a 2-dimensional  cellular  array,  was  designed  to  exhibit 
one  of  the  major  features  of  life:  auto-reproduction.  As  a mathematician,  von  Neumann 
perceived  the  theory  of  cellular  automata  as  a mathematical  tool  to  study  parallel  machines: 
“understanding  their  mathematical  structure  will  prepare  us  eventually  to  go  further  into 
the  almost  unexplored  theory  of  parallel  computers”  (page  xi)  [8]. 
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J.  Von  Neumann’s  work  on  cellular  automata  has  inspired  a series  of  very  interesting 
research  in  computer  architecture,  computation  theory,  computer  modeling,  and  applica- 
tions [9,10,11,12,13,14].  Among  others,  Unger  was  the  first  to  suggest  the  use  of  a cellular 
array  machine  for  parallel  image  processing  [15,16].  Since  then,  many  special  parallel  image 
processing  computers,  such  as  the  CLIP  series  [17,18,19]  and  the  Cytocomputer  [20,21], 
have  been  developed  and  have  been  shown  to  be  very  successful. 

Along  with  the  development  of  the  cellular  architectures,  some  special  algebraic  systems 
for  image  processing  have  been  proposed.  G.  Matheron  and  J.  Serra  formulated  an  alge- 
braic theory  known  as  mathematical  morphology  [22,23,24],  which  was  also  independently 
investigated  by  Sternberg  [25,26].  Mathematical  morphology  is  about  the  study  of  shapes 
and  patterns  in  general  by  the  use  of  structure  elements.  It  is  based  on  the  operations 
of  Minkowski  addition  and  subtraction  of  sets  in  n-dimensional  Euclidean  space  Rn.  The 
Minkowski  operations  are  often  referred  to  as  the  morphological  dilation  and  erosion  oper- 
ations. It  turns  out  that  these  morphological  operations  and  techniques  can  be  applied  to 
many  image  processing  and  image  analysis  problems. 

The  morphological  approach,  however,  has  its  limitations  in  expressing  some  commonly 
used  image  processing  algorithms  such  as  linear  image  to  image  transformations  [27].  To 
remedy  the  limitations  of  the  morphological  approach  and  establish  a universal  image  pro- 
cessing system,  G.  Ritter  at  the  University  of  Florida  started  the  development  of  the  image 
algebra,  which  has  been  accepted  as  the  Air  Force  Standard  Image  Algebra  [1].  It  has  been 
shown  that  the  image  algebra  is  capable  of  expressing  all  image  to  image  transformations 
[2].  In  particular,  it  incorporates  linear  algebra  and  minimax  algebra  as  subalgebras  [28,29]; 
therefore  it  can  be  used  to  express  all  linear  and  morphological  transformations. 
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1.2  Recursive  Methods  in  Signal  and  Tmage  Processing 

Many  useful  image  transformations,  such  as  the  Fourier  transform,  average  and  median 
filters,  Kirsch  and  some  other  edge  detectors,  are  considered  to  be  parallel  image  operations. 
Each  parallel  image  operation  can  be  performed  independently  on  each  pixel  of  a given  image 
since  each  pixel  value  of  the  transformed  image  is  only  a function  of  some  pixel  values  of 
the  given  image.  The  sequence  in  which  the  pixels  are  processed  is  completely  irrelevant; 
therefore,  these  parallel  image  operations  can  be  applied  to  each  pixel  simultaneously  if 
appropriate  parallel  hardware  is  available.  Parallel  image  transformations  are  also  referred 
to  as  nonrecursive  transformations. 

Along  with  nonrecursive  transformations,  a class  of  recursive  transformations  are  also 
widely  used  in  signal  and  image  processing,  e.g.,  HR  (infinite  impulse  response)  filters, 
sequential  block  labeling,  AR  (auto- regressive)  image  modeling,  predictive  coding,  adaptive 
dithering,  etc.  [30,31,32,33,34,35].  One  of  the  characteristic  properties  of  these  recursive 
image  operations  is  that  a pixel  value  of  the  transformed  image  may  depend  on  both  the 
pixel  values  of  the  given  image  and  some  of  the  pixel  values  of  the  transformed  image 
itself.  Thus,  the  transformed  image  may  have  to  be  computed  recursively  according  to 
some  partial  order  imposed  on  the  underlying  image  domain.  In  other  words,  a pixel  value 
of  the  transformed  image  may  not  be  processed  until  all  the  pixels  ordered  previous  to 
it  have  been  processed.  Some  of  these  recursive  operations  have  a parallel  counterpart, 
and  some  of  them  are  inherently  sequential  in  nature.  Yet,  many  recursive  operations  are 
considered  to  be  more  efficient  if  only  conventional  sequential  computers  are  employed  [36]. 

A study  of  the  recursive  image  operations  in  the  context  of  a coherent  algebraic  system 
such  as  the  image  algebra  is  not  only  desirable  but  also  necessary  to  exhibit  the  expressive 
power  and  usefulness  of  the  image  algebra.  In  order  to  specify  these  recursive  transfor- 
mations in  the  image  algebra,  some  extra  operations  must  be  defined  and  their  algebraic 
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properties  investigated  within  the  framework  of  the  image  algebra  so  that  a better  under- 
standing and  implementation  of  recursive  image  operations  may  be  achieved. 

1.3  Summary  of  Results 

The  main  result  of  this  research  is  a theoretical  framework  based  on  the  image  algebra 
for  studying  recursive  operations  in  image  processing.  First,  the  notions  of  recursive  tem- 
plates and  recursive  template  operations  are  introduced  and  computation  schemes  for  these 
operations  are  discussed.  The  recursive  template  operations  are  capable  of  representing 
linear  and  nonlinear  recursive  transformations  commonly  encountered  in  image  processing. 
Their  algebraic  properties  are  then  stated  and  proved,  providing  a mathematical  basis  for 
recursive  template  composition  and  decomposition.  A major  contribution  of  this  research 
is  the  development  of  the  relationships  between  the  extended  image  algebra  and  other  al- 
gebraic structures.  In  the  linear  case,  the  extended  image  algebra  may  be  mapped  to  rings; 
in  the  nonlinear  case,  it  may  be  mapped  to  a subalgebra  of  the  minimax  algebra.  Another 
group  of  results  concerns  the  problem  of  template  decomposition.  Template  decomposition 
plays  an  important  role  in  image  algorithm  optimization  and  parallel  image  processing. 
Some  necessary  and  sufficient  conditions  along  with  algorithms  are  given  for  decomposing 
recursive  and  nonrecursive  templates. 

The  organization  of  the  dissertation  is  as  follows: 

An  overview  of  the  image  algebra  and  related  subjects  such  as  recursive  filters,  the  z- 
transform,  and  the  minimax  algebra  are  given  in  Chapter  2.  Basic  definitions  and  results 
from  previous  relevant  research  are  given  and  will  be  referenced  in  the  rest  of  the  disserta- 
tion. In  Chapter  3,  we  introduce  the  new  definitions  of  recursive  templates  and  operations 
into  the  current  body  of  the  image  algebra.  These  recursive  template  operations  are  de- 
fined as  image  to  template  operations  as  well  as  template  to  template  operations.  Different 
computation  schemes  for  the  recursive  operations  will  also  be  discussed.  In  Chapter  4,  we 
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study  the  algebraic  properties  of  these  recursive  operations  and  their  relationships  to  other 
algebraic  structures.  We  also  introduce  the  max-min  transforms  and  show  that  they  are 
useful  to  describe  recursive  morphological  filters.  In  Chapter  5,  we  give  some  necessary  and 
sufficient  conditions  for  decomposing  certain  classes  of  recursive  and  nonrecursive  templates. 
Algorithms  for  decomposing  these  templates  will  be  described  in  detail.  In  Chapter  6,  we 
present  applications  of  the  recursive  template  operations  to  signal  and  image  processing. 
Finally,  we  give  concluding  remarks  and  suggestions  for  further  research. 


CHAPTER  2 
BACKGROUND 


This  chapter  contains  the  basic  concepts  and  definitions  of  the  image  algebra.  It  also 
gives  a brief  review  of  recursive  filters  and  z-transforms  as  well  as  the  minimax  algebra. 
The  notation  presented  here  will  be  used  throughout  this  dissertation. 

2.1  Image  Algebra 

The  image  algebra  is  an  abstract  algebraic  structure  specially  designed  for  image  pro- 
cessing [1].  In  this  structure,  basic  objects  in  image  processing  such  as  images,  point  sets, 
and  templates  are  formally  defined  as  its  operands.  A set  of  unary  and  binary  operations 
are  defined  on  and  between  these  operands  to  express  various  image  processing  transforms. 
It  has  been  demonstrated  that  many  commonly  used  image  processing  transformations, 
such  as  generalized  convolutions,  DFT,  edge  detectors,  and  morphological  operations,  can 
be  easily  expressed  within  the  framework  of  the  image  algebra  [2,3,28,4,29], 

For  the  purpose  of  this  dissertation,  only  basic  definitions  and  operations  of  the  image 
algebra  are  given.  A complete  description  of  the  image  algebra  can  be  found  in  Ritter  et 
al.  [1]. 

2.1.1  Images 

Images  are  a basic  type  of  operand  in  the  image  algebra.  Images  are  defined  in  terms 
of  two  other  elementary  types  of  operands,  namely,  value  sets  and  point  sets. 

A value  set  is  a homogeneous  algebra.  The  most  commonly  used  value  sets  in  image 
processing  are  the  set  of  integers,  real  numbers,  extended  real  numbers  (which  include  one 
or  both  of  the  symbols  — oo  or  +oo),  complex  numbers,  and  binary  numbers  of  fixed  length 
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k,  which  will  be  denoted  by  Z,  R,  R_00,  R+0O,  C,  and  Z 2k,  respectively.  We  denote  an 
arbitrary  value  set  by  F. 

A point  set  is  a subset  of  a topological  space.  Point  sets  are  usually  denoted  in  the 
image  algebra  by  the  letters  X,  Y and  W,  and  their  elements  by  bold  lower  case  letters, 
x,y,  etc..  Rectangular  and  hexagonal  discrete  arrays  are  usually  the  assumed  point  sets  in 
digital  image  processing. 

Given  point  and  value  sets  X and  F,  respectively,  an  F valued  image  a on  X is  the 
graph  of  a function  a : X — ► F.  Thus,  an  F valued  image  a on  X is  of  the  form: 

a = {(x,a(x)) : x G X} 

where  a(x)  G F. 

The  set  X is  called  the  set  of  image  points  of  a,  and  the  range  of  the  function  a (which 
is  a subset  of  F)  is  the  set  of  image  values  of  a.  The  set  of  all  F valued  images  on  X is 
denoted  by  Fx. 

2.1,2  Generalized  Templates 

Another  basic  type  of  operand  in  the  image  algebra  is  generalized  template,  which  plays 
a very  important  role  in  expressing  all  kinds  of  image  transformations. 

Let  X and  Y be  point  sets  and  F a value  set.  A generalized  F valued  template  t from  Y 
to  X is  a function  t : Y — ► Fx.  Thus,  for  each  y G Y,  t(y)  is  an  F valued  image  on  X.  For 
notational  convenience  we  define  ty  = t(y).  The  sets  Y and  X are  called  the  target  domain 
and  range  space  of  t,  respectively.  The  point  y is  called  the  target  point  oft  and  the  values 
ty(x)  are  called  the  weights  of  the  template  t at  y.  The  set  of  all  F valued  templates  from 
Y to  X will  be  denoted  by  (FX)Y. 

If  t is  a real  valued  template  from  Y to  X,  then  we  define  the  support  of  ty  to  be 
S(ty)  = {x  G X : ty(x)  # 0}.  If  t is  an  extended  real  valued  template,  then  we  define 
5-oo(ty)  = {x  G X : ty(x)  / -oo}  and  S+00(ty)  = {x  G X : ty(x)  7*  +00}. 
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A generalized  template  in  image  algebra  is  a generalization  of  the  concepts  of  template, 
mask,  window,  and  structuring  element  as  used  in  image  processing. 

An  important  property  of  the  generalized  templates  is  their  translation  invariance.  A 
template  t £ (Fx)x  is  called  translation  invariant  if  and  only  if  for  each  triple  x,y  £ X 
with  y + z £ X and  x + z £ X,  we  have  ty(x)  = ty+z(x  + z).  Image  transformations  with 
invariant  templates  may  be  easily  mapped  to  SIMD  and  systolic  parallel  architectures. 

To  facilitate  the  descriptions  of  some  image  transformations,  the  transpose  and  dual  of 
certain  templates  are  defined  in  the  image  algebra.  The  transpose  of  a template  t £ (FX)Y 
is  the  template  t'  6 (FY)X  defined  by  t^(y)  = ty(x).  If  t £ (RX00)Y  (or  t £ (R?oo)Y), 
then  the  additive  dual  of  t is  the  template  t*  £ (R^co)x  (or  t*  £ (R^^)  defined  by 
tx(y)  = — ty(x).  If  t G ((R^)X)Y,  then  the  multiplicative  dual  of  t is  the  template 
t € ((R-$^)Y)X  defined  by 


tx(y) = < 


l/ty(X) 

+oo 

0 


if  ty(x)  ^ 0 and  ty(x)  ^ +oo 
if  ty(x)  = 0 
if  ty(x)  = +oo 


2.1.3  Operations  on  Images 


The  basic  operations  on  F valued  images  are  derived  from  the  naturally  induced  oper- 
ations on  the  value  set  F.  In  particular,  the  basic  binary  operations  of  addition,  multipli- 
cation, and  maximum  on  Rx  are  defined  as  follows: 

If  a,  b G RX,  then 


a + b = {(x,  c(x)) : c(x)  = a(x)  + b(x),x  G X} 

a*b  = {(x, c(x)) : c(x)  = a(x)  • b(x), x £ X} 

a V b = {(x,  c(x)) : c(x)  = a(x)  V b(x),  x G X}. 

The  exponentiation  of  two  real  images  a,  b G RX  is  defined  as 


ab  = {(x,  c(x))  : c(x)  = a(x)b^xi  if  a(x)  ± 0,  otherwise  c(x)  = 0,  x G X} 
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If  X is  finite,  then  the  dot  product  of  two  images  a,  b G RX  is  defined  as 


a • b = ^2  a(x)  ' b(x) 
x€X 

The  additive  dual  of  an  extended  real  valued  image  a G R*^  ( or  a G R^J,  denoted 
by  a*,  is  defined  by  a*(x)  = -a(x),  where  we  define  — (-00)  = +00  and  —(+00)  = -00. 

If  a G (R^)X,  where  R^  = {r  G R : r > 0}  U {+00},  then  the  multiplicative  dual  of 
a,  denoted  by  a,  is  defined  by 


a(x) 


l/a(x)  if  a(x)  7^  0 and  a(x)  ^ +00 
< +00  if  a(x)  = 0 

0 if  a(x)  = +00 


Characteristic  functions  are  useful  operations  in  the  image  algebra.  If  a,  b G RX,  then 
the  characteristic  function  x<b(a)  is  defined  by 


X<b(a)  = {(x,c(x))  : c(x)  = 1 if  a(x)  > b(x),  else  c(x)  = 0}. 

Other  characteristic  functions  x<b>  X=b>  X>b>  and  X>b>  are  defined  similarly. 

2.1.4  Operations  between  Images  and  Templates 

Operations  between  images  and  templates  are  considered  to  be  higher  level  operations 
than  the  operations  between  images  mentioned  above  since  most  global  and  neighborhood 
image  transformations,  such  as  convolution,  DFT,  feature  extraction,  and  morphological 
transformations,  can  be  easily  described  with  these  image-template  operations. 

First  a global  reduce  function  T : Fx  -*  F is  defined  as  follows.  Let  X C R"  be  finite, 
X = {x1?x2, . • .,xm}.  If  7 be  an  associative  and  commutative  binary  operation  on  the 
value  set  F , then  the  global  reduce  operation  T on  Fx  induced  by  7 is  defined  by 


Ta  = rxeXa(x)  = a(xi)7a(x2)7  • • -7a(xm) 


where  a G Fx. 
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Now,  let  Fj  and  F2  be  two  value  sets  and  o : Fj  x F2  ->F  a binary  operation.  If  7 
is  an  associative  and  commutative  binary  operation  on  F,  a 6 Fx,  and  t 6 (FX)Y,  then 
the  generalized  right  product  of  a with  t is  the  binary  operation  © : Fx  x (FX)Y  — ► FY 
defined  by 

a©t  = {(y,b(y)):b(y)  = rxeXa(x)oty(x),y  G Y}. 

Ift  G (FY)X,  0':F2xF1-4  F,  then  the  generalized  left  product  ©'  : (FY )X  X Fx  — > FY 
is  defined  by 

t ©'a  = {(y,  b(y))  : b(y)  = rxeXtx(y)  o'  a(x),  y G Y}. 

Although  one  can  derive  many  image-template  operations  by  substituting  appropriate 
binary  operations  for  7 and  o in  the  above  definitions  of  the  generalized  forward  and  back- 
ward template  operations,  there  are  three  basic  image-template  operations  in  the  image 
algebra,  which  are  sufficient  to  describe  most  commonly  used  image  transforms.  The  three 
basic  operations  are  denoted  by  ©,0,  and  ©,  and  are  called  generalized  convolution,  ad- 
ditive maximum,  and  multiplicative  maximum,  respectively.  The  operation  © is  defined  for 
real  and  complex  valued  images  while  the  operations  0 and  © may  be  applied  to  extended 
real  valued  images. 

Let  X C R"  and  Y C Rm.  Suppose  a G RX,  t G (RX)Y.  Then  the  generalized  convo- 
lution is  defined  as 


a©t  = {(y,  c(y))  : c(y)  = ^ a(x)  • ty(x),  y G Y}. 

xeX 

If  a G RXoo  and  t G (RXQO)Y,  then  the  additive  maximum  is  defined  as 


a®t  = {(y,c(y)):c(y)=  V a(x)  + ty(x),  y G Y}, 

xeX 

where  a(x)  + ty(x)  = max{a(x)  + ty(x)  : x G X}.  We  assume  the  extended  arithmetic 
xeX 

convention  that  r + (-00)  = (-00)  + r=  -00  for  all  r G R-oo- 
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Similarly,  if  a G (Rp^)X  and  t E ((R^)X)Y,  then  the  multiplicative  maximum  is 
defined  as 

a©t  = {(y,c(y)):  c(y)  = \J  a(x)  • ty(x),  y G Y}, 

xeX 

where  we  assume  the  extended  arithmetic  convention  that  r • (+oo)  = (+oo)  • r = -foo  for 
all  r > 0,  and  0 • (+<x>)  = (+oo)  -0  = 0. 

The  dual  operations,  additive  minimum,  denoted  by  O,  and  multiplicative  minimum, 
denoted  by  ©,  can  be  defined  in  terms  of  the  basic  template  operations  and  the  dual  of 
images  and  templates  as  follows: 

aEJt  = (t*Ma*)* 

and 

a@t  = (t©a). 

2.1.5  Operations  between  Generalized  Templates 

The  basic  operations  between  generalized  templates  are  defined  in  terms  of  the  basic 
operations  on  images.  Let  s,t  G (RX)Y.  The  addition,  multiplication,  and  maximum  of 
two  templates  are  defined  as  follows: 


s + t 

= {(y>ry): 

: ry  = sy  + ty,  y e Y) 

s * t 

= {(y,ry): 

ry  = sy  ‘ ty>  y £ Y} 

s V t 

= {(y>ry): 

: ry  = sy  v , y G Y} 

Similarly,  the  high  level  operations  ©,0,  and  © between  images  and  templates  gen- 
eralize to  operations  between  templates.  The  use  of  these  high  level  operations  to  combine 
two  templates  is  often  referred  to  as  template  composition.  Let  t E (RX)Y  and  s E (Rw)x 
be  two  real  valued  templates.  The  template  composition  r = s©t  is  a new  template  in 
(Rw)y,  and  is  defined  by 


r — {(y>  ry)  : ry  — ty  ©s',  y G Y}, 
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where  ry  = ty©s'  is  the  regular  generalized  convolution  between  image  and  template.  In 
other  words, 

ry  = {(w>ry(w))  : ry(w)  = ^ ty(x)-s^(x),  w £ W} 

xeX 

= {(w,  ry(w))  : ry(w)  = ^ ty(x)  • sx(w),  w £ W}. 

xeX 

Note  that  if  t and  s are  invariant  templates,  then  r would  be  an  invariant  template  too. 
In  this  case,  to  get  the  new  template  r,  one  needs  only  to  compute  ry  once  for  any  y £ Y. 

The  operations  0 and  © between  templates  are  defined  in  terms  of  the  corresponding 
image-template  operations  also.  If  t G (R?:00)Y  and  s € (R^c)x,  then  r = s [V] t £ (R^Q)Y 
is  defined  by 

r = {(y>ry)  : ry  = tyE]s',  y £ Y}. 

If  t £ ((Rp^)x)Y  and  s £ ((R^)w)x,  then  r = s©t  £ ((Rp^)W)Y  is  defined  by 

r = i(y>  ry)  : ry  = ty  ©s',  y £ Y}. 

The  dual  operations  E3  and  © between  templates  are  given  by 

sElt  = (t*(V]s*)* 

and 

s@t  = (t©s). 

Many  other  image  and  template  operations  are  described  in  Ritter  et  al.  [1],  Some 
extensive  mathematical  studies  of  these  template  operations  can  be  found  in  Gader  [28] 
and  Davidson  [29]. 

2.2  Recursive  Filters  and  z-transform 

Two-dimensional  digital  filters  have  been  studied  and  applied  to  image  processing  ap- 
plications such  as  image  enhancement  and  restoration  for  some  time  [31,37,30,35].  In  this 
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section,  The  definitions  and  basic  properties  of  recursive  filter  and  the  z-transform  are 
presented. 

2.2.1  Two-dimensional  Signals  and  Systems 

A two-dimensional  (2-D)  discrete  space  signal  or  sequence  can  be  represented  by  a 
function  of  two  integer-type  arguments.  The  2-D  unit  impulse  sequence,  denoted  by  <5(m,  n), 


is  defined  as 


6(m,  n ) — 


1,  if  m — n — 0 


0,  otherwise. 

A 2-D  sequence  x(m , n)  is  called  a separable  sequence  if  it  can  be  expressed  in  the  form 

x(m,  n)  = /(m)  ■ g(n), 

where  /(m)  and  g(n ) are  single- variable  functions.  A 2-D  sequence  x(m,n ) is  called  causal 
if  x(m,  n)  = 0 for  m < 0 or  n < 0. 

A 2-D  discrete  system  or  transformation  L is  described  by  the  equation 

y(m,  n ) = L[x(m , n)] 

where  x(m,n)  and  y(m,n ) are  2-D  input  and  output  sequences,  respectively. 

A 2-D  system  L is  linear  shift-invariant  if  it  satisfies  the  following  relations: 

L[a\X\{m,  n)  + 02X2(m,  n)]  = a\L[x\(m,  n)]  + a,2L[x2(m , n)] 

L[x(m  — mo,  n — no)]  = y(m  — mo,  n — no) 

A 2-D  linear  shift-invariant  (LSI)  system  may  be  characterized  by  its  unit  impulse 
response  h(m,n)  in  terms  of  2-D  convolution  [31] 


oo  oo 


y(m,n ) = L[x(m,n)\  = x(m,n)®h(m,n ) = E E h(k,l)x(m  — k,n  — l), 

k=— oo  (=— oo 

where  © is  the  convolution  operator  and  fi(m,  n)  = L[S(m,  n)],  the  system  response  of  the 
unit  impulse  sequence. 
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A 2-D  linear  shift-invariant  system  is  called  separable  if  its  unit  impulse  response  h(m,  n ) 
is  a separable  sequence.  Separable  systems  play  an  important  role  in  analysis  and  compu- 
tation of  2-D  systems  as  we  will  see  later. 

The  system  transfer  function  H(zx,z2 ) is  defined  by  the  z-transform  of  h(m,n ) 

oo  oo 

H{z1,z2)  = H h(m,n)zimZ2n 

m~ — oo  n=— oo 

which  is  equivalent  to  the  unit  impulse  response  h(m,  n)  in  defining  an  LSI  system. 

2.2.2  z- Transform 

The  z-transform  can  be  used  to  represent  a large  class  of  sequences,  including  many 
unstable  ones.  The  z-transform  method  provides  a useful  technique  for  transforming  a 
difference  equation  into  an  algebraic  equation  so  that  well-defined  mathematical  analysis 
may  be  applied  to  the  underlying  system. 

The  z-transform  of  a sequence  x(m,n ) is  denoted  by  X(zx,z2)  and  is  defined  by 

OO  OO 

X(zi,z2)=  ^ x(m,n)z^mz^n 

m=— oo  n=— oo 

where  zx  and  z2  are  complex  variables. 

Let  Xi(z\,z2 ) and  X2(zi,z2)  be  the  z-transform  of  the  input  sequences  xx (m,n)  and 
x2(m,n),  respectively.  Some  useful  properties  of  the  z-transform  are: 

1.  Linearity  property 

axa:i(m,  n)  + a2x2(m,n ) o axXi(zi,z2)  + a2X2(zx,z2) 

2.  Shifting  property 

mo  — 1 no  — 1 

xx{m  — m0,  n — n0)  zxm°  z^710  X\(zx,  z2)  + xx(k  - m0,  / - nQ)zxkz^1 

k= o l=o 

3.  Convolution  property 


x1(m,n)®x2(m,n)  = Xx(zi,  z2)  ■ X2(zu  z2) 
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4.  Separability  property 


*1  («*)  • *2(»)  = ^i(^i)  • AT2(z2) 


By  the  convolution  property  of  the  z-transform, 

which  is  another  way  to  define  the  system  transfer  function  H(zi,  z2). 

2.2.3  Recursive  Filters 

A 2-D  recursive  filter  is  characterized  by  a 2-D  difference  equation 

y(m,n)=  ^ a(k,l)x(m  - k,n  - l)  - ^ b(k,l)y(m  - k,n  - l) 

k,l€Ra  k,leRb-(  0,0) 

where  x(m,n)  and  y(m,n ) are  input  and  output  sequences,  Ra  and  Rb  are  the  regions  of 
support  for  constants  a(k,l ) and  b(k,l),  respectively. 

Applying  the  z-transform  to  the  difference  equation  and  using  the  shifting  property, 

Y(zuz2)=  J2  a(k,l)zikzJlX(z1,z2)-  Kk^)zikz2lY{zuz2) 

k,ieRa  k,ieRb-(o,o) 

So  the  system  transfer  function  of  a 2-D  recursive  filter  is  of  the  form 

H(zuz2)  = Y(Zl’Z2)  = a(fc»  0*i~V 

X (2l’  ~2)  1 + ^2k,ieRb-(o,o)  l)zx  kz2  1 

A 2-D  recursive  filter  is  called  quarter-plane  ( QP ) if  Rb  is  in  one  of  the  four  quarter 
planes. 

If  the  region  of  support  Rb  is  empty,  then 

a^)  = Ytr^=  E <«u)*rV 

A(~i  ,z2)  ki€Ra 

which  defines  a nonrecursive  filter. 

A filter  is  called  a simplified  recursive  filter  if  its  system  transfer  function  is  of  the  form: 


#(*!,**)=  = 


X(21’  1 + ^2k,leRb-(0, 0)  021  kz2  1 
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Note  that  the  system  transfer  function  of  a recursive  filter  can  be  expressed  as 


h(z1,Z2)=  = 


--  a(k,l)z1kz2l 


X(zi,Z2)  1 + J2k,ieRb-(o,o)  0zi  kz2  1 k,ieRa 
which  is  the  cascade  of  a simplified  recursive  filter  and  a nonrecursive  filter.  This  function 
is  also  referred  to  as  the  canonical  direct  form  realization  [31]. 


2.3  Ordered  Algebraic  Structures 

Ordered  algebraic  structures  such  as  lattice,  lattice-ordered  groups,  etc.,  have  been 
developed  in  mathematics  to  deal  with  many  operation  research  and  combinatorial  opti- 
mization problems.  Cuninghame- Green  systematically  presented  such  an  algebraic  system 
called  minimax  algebra  [38],  which  include  structures  on  matrices.  In  this  section,  we  briefly 
review  some  related  concepts  and  properties  of  the  minimax  algebra. 

Let  F be  a linearly  ordered  set  and  V a binary  operation  on  F.  The  structure  (F,  V)  is 
called  a commutative  band  or  semi-lattice  if  it  satisfies  the  following  axioms: 


x V (y  V z)  = (x  V y)  V z 

(2.1) 

x V y — y V x 

(2.2) 

x V x = x 

(2.3) 

for  all  x,y,z  6 F. 

Let  X be  a binary  operation  on  F.  The  structure  (F,  V,  x)  is  called  a belt  if  it  satisfies 
the  following  properties  in  addition  to  the  three  properties  above: 


x x (y  x z)  = (x  x y)  x z 

(2.4) 

x X (y  V z)  = (x  X y)  V {x  X z) 

(2.5) 

(y  V 2)  x x = (y  x x)  V (z  x x) 


(2.6) 
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for  all  x,y,z  £ F.  The  operation  V is  called  an  addition,  and  the  operation  X is  called  a 
multiplication. 

An  element  0 £ F is  called  zero  or  null  element  if 

x V 0 = x 

xx0=0xx=0 

for  all  x £ F. 

An  element  1 £ F is  called  identity  or  unit  element  if 

xxl^lxx  — x 

for  all  x £ F. 

If  (F,  x)  is  a commutative  semigroup,  then  the  structure  (F,  V,  x)  is  called  a commu- 
tative belt.  If  (F  - {0},x)  is  a group,  then  the  structure  (F,V,x)  is  called  a division 
belt. 

Now,  let  A and  x'  be  two  more  binary  operations  on  F.  The  structure  (F,  V,  x,  A,  x')  is 
called  a belt  with  duality  if  both  (F , V,  x ) and  (F,  A,  X7)  are  belts  and  the  following  property 
holds 


x V (y  A a:)  = x A (y  V x)  = x (2.7) 

for  all  £,?/,  z £ F.  The  operation  A is  called  a dual  addition  and  the  operation  x'  is  called 
a dual  multiplication. 

Some  examples  of  belt  with  duality  are  (R_oo,  V, +,  A, +)  and  (R-°,  V,  *,  A,  *). 

Let  Fi  and  F2  be  two  belts.  A function  / : Fj  — ► F2  is  called  a belt  homomorphism  if 


it  satisfies 


f(x  Vy)  = f(x)  V f(y) 
f(x  xy)  = f(x)  x f(y) 
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for  all  x,y  £ Fi.  If  the  function  / is  1-1  and  onto,  then  it  will  be  an  isomorphism. 

The  belt  structures  may  be  defined  on  matrices  also.  Let  (F,  V,  x)  be  a belt  and  Mn  be 
the  set  of  all  n X n matrices  with  values  in  F.  Let  A = (a,y),  B = (6,y)  £ Mn.  We  define  a 
matrix  addition  V and  a matrix  multiplication  x on  Mn  as  follows: 


(aij)  V (bij)  — (dij  V bij ) 

n 

( aij ) X (Pi j)  = ( \J  aik  X bkj ) 

A:=l 

Similarly,  we  can  define  a dual  matrix  addition  A and  a dual  matrix  multiplication  x' 


by 


( d ij ) A (b^)  — (dij  A bij) 


( aij ) X (bij)  — ( \J  dik  X bkj) 
k= l 


It  has  been  shown  [38]  that  the  structure  (Mn,  V,  x)  is  a belt  and  (Mn,  V,  x,  A,  x')  is  a 
belt  with  duality. 

In  this  research,  we  are  interested  in  the  following  two  matrix  belts 
• the  matrix  belt  (Mn,  V,  +)  derived  from  the  belt  (R_oo,  V,  +),  and 


• the  matrix  belt  (Mn,  V,  *)  derived  from  the  belt  (R^°,  V,  *) 

Note  that  -oo  and  0 are  the  zero  and  identity  elements  of  the  belt  (R^,  V,  +),  respec- 
tively; 0 and  1 are  the  zero  and  identity  elements  of  the  belt  (R-°,  V,*),  respectively. 
Thus,  for  the  matrix  belt  (Mn,  V,  +),  we  have  the  zero  matrix 

— OO  . . . —00 

— oo  . . . — oo 

— oo  . . . — oo 

— oo  ...  — oo 
0 ...  — oo 

—oo  ...  0 


— oo 
— oo 


and  the  identity  matrix 


On  = 


— oo 
0 

— oo 


-oo 
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For  the  matrix  belt  ( Mn , V,  *),  we  have  the  zero  matrix 


and  the  identity  matrix 


0n 


0 0 
0 0 


0 

0 


0 0 ...  0 


1 0 ...  0 

0 1 ...  0 

0 0 ...  1 


CHAPTER  3 

RECURSIVE  TEMPLATES  AND  OPERATIONS 


In  this  chapter,  we  introduce  the  notion  of  generalized  recursive  templates  and  recursive 
template  operations,  which  are  the  direct  extension  of  the  generalized  templates  and  the 
corresponding  template  operations  defined  in  the  image  algebra.  First,  we  begin  with  the 
definition  of  recursive  templates.  Next,  we  introduce  the  basic  recursive  image-to-template 
operations  and  discuss  the  possible  parallel  computation  of  these  operations. 

3.1  Generalized  Recursive  Templates 

A generalized  recursive  template  is  defined  in  terms  of  the  regular  generalized  template 
given  in  Chapter  2 and  some  partial  order  imposed  on  the  target  point  set  involved. 
Definition  A partially  ordered  set  {P,<)  (or  poset)  is  a set  P together  with  a binary 
relation  X,  satisfying  the  following  three  axioms: 

1.  For  all  x 6 P,  x -<  x.  ( reflexivity ) 

2.  If  x -<  y and  y -<  x,  then  x = y.  ( antisymmetry ) 

3.  If  x -<  y and  y -<  z,  then  x -<  z.  ( transitivity ) 

Definition  Let  (P, -<)  be  a poset,  and  [n]  = {1,2, where  n - card(P).  A linear 
extension  of  (P,  -<)  is  a bijection  a : P ->  [n]  such  that  x -<  y in  P if  and  only  if  a(x)  < a{y) 
in  [rc]. 

Let  X and  Y be  point  sets  and  F a value  set.  Also,  let  -<  be  a partial  order  imposed 
on  the  point  set  Y.  A generalized  F valued  recursive  template  t from  Y to  X is  a function 
t : Y -4-  (Fx,Fy),  where  t/  : Y Fx,  and  tx  : Y — ► FY,  such  that 
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1-  y £ <S(My))  and 

2.  for  each  z £ «S(t^(y)),  z -<  y. 

That  is,  the  Support  of  is  consistent  with  the  partial  order  imposed  on  Y.  Thus,  for 
each  y £ Y,  t^(y)  is  an  F valued  image  on  X,  and  t^(y)  is  an  F valued  image  on  Y. 

Note  that  the  consistency  requirement  for  is  important  because  it  ensures  the  recur- 
sive computability  of  the  defined  template.  In  fact,  a recursive  template  could  be  consistent 
with  many  different  partial  orders  on  Y,  and  any  of  their  linear  extensions.  Hereafter,  when 
we  define  a recursive  template  that  is  consistent  with  some  partial  order  -<,  we  will  say  that 
the  recursive  template  admits  the  partial  order  -<. 

For  our  convenience,  we  define  t*y  = t/(y),  and  txy  = tx(y).  So,  ty  = (t^y,t^y). 

As  the  case  of  the  regular  generalized  templates,  the  sets  Y and  X are  called  the  target 
domain  and  range  space  of  t,  respectively.  The  point  y is  called  the  target  point  oft,  and  the 
values  t^y(x)  and  t^y(x)  are  called  the  weights  of  recursive  part  and  nonrecursive  part  of 
the  template  t,  respectively.  The  set  of  all  F valued  recursive  templates  from  Y to  X that 
admit  the  partial  order  on  Y will  be  denoted  by  (FX,F^)^,  or  simply  by  (Fx,F^r)^r. 

The  support  of  a recursive  template  is  defined  similarly  to  the  one  for  a regular  template. 
In  particular,  if  t is  a real  valued  recursive  template  from  Y to  X,  then  we  define  the  support 
of  ty  to  be  <S(ty)  = (5(t^y),<S(t^y)).  If  t is  an  extended  real  valued  recursive  template, 
then  we  define  S-co(ty)  = (S.^t^S.^y)),  and  S+0O(ty)  = («S+00(t^y),<S+00(t^y)). 
Note  that  if  the  support  S(t^y)  is  empty  for  each  y £ Y,  then  t would  be  collapsed  to  a 
nonrecursive  generalized  template. 

A recursive  template  t £ (FX,FX)X  is  called  translation  invariant  if  and  only  if  for 
each  triple  x,  y £ X with  y + z £ X and  x + z £ X,  we  have  ty (x)  = ty_|_z(x  + z),  or  equiv- 
alently, t/y(x)  = t^y+z(x  + z)  and  txy(x)  = txy+z(x  + z). 

An  example  of  an  invariant  recursive  template  t is  given  below: 
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t-t  = 


If  t is  an  invariant  recursive  template  and  it  has  only  one  pixel  value  defined  on  the  target 
point  of  its  nonrecursive  support  S(i^),  than  t is  called  a simplified  recursive  template. 
Pictorially,  a simplified  recursive  template  can  be  drawn  in  the  same  way  as  a nonrecursive 
template  since  the  recursive  part  and  the  nonrecursive  part  do  not  overlap.  So,  we  may 
draw  the  above  template  t,  which  is  a simplified  recursive  template,  as  follows: 


-1 

1 

1 

s N 
w 

The  transpose  and  dual  for  recursive  templates  are  defined  in  terms  of  those  for  non- 
recursive  templates.  The  transpose  of  a recursive  template  t £ (FX,FX)X  is  the  template 
t'  £ (FX,FX)X  defined  by  t'  = (t^,t^).  Similarly,  if  t 6 (R^,  R?00)X,  then  the  additive 
dual  of  t is  defined  by  t*  = (tj,t^).  If  t £ ((R|^)x,  (R$°  )X)X,  then  the  multiplicative 
dual  of  t is  defined  by  t = (t^,t_-). 

3.2  Operations  between  Images  and  Recursive  Templates 

There  are  three  basic  template  operations  defined  in  the  image  algebra  that  are  used  to 
express  nonrecursive  transformations  of  real  valued  (and/or  extended  real  valued)  images. 
They  are  generalized  convolution,  additive  maximum , and  multiplicative  maximum , denoted 
by  ©,0,  and  ©.  In  order  to  express  recursive  transforms  in  the  image  algebra,  we  define 
three  new  recursive  template  operations,  denoted  by  ®^,  0^,,  and  @x,  respectively.  They 
operate  between  images  and  recursive  templates. 


= 


-l 

l 

l 

V J 
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Let  Y,X  C Rn  be  finite.  Let  a E Rx  and  t E (RX,RY  )Y.  We  define  the  generalized 
recursive  convolution  as 

a®xt  = {(y,c(y))  : c(y)  = ^ a(x)  • t/y(x)  + c(z)  • t-<y(*)»  y € Y}. 

x€<S(t,<y)  Z65(t_<y) 

The  recursive  template  operations  compute  a new  pixel  value  c(y)  based  on  both  the 
pixel  values  a(x)  of  the  given  image  and  some  previously  calculated  new  pixel  values  c(z), 
which  are  determined  by  the  given  partial  order  -<  and  the  region  of  support  of  the  partici- 
pating template.  By  the  definition  of  the  recursive  templates,  z -<  y for  all  z £ <S(t^y)  and 
y £ <S(txy).  Thus,  c(y)  is  always  recursively  computable.  Some  commonly  used  partial  or- 
ders in  2-D  recursive  transforms  are  forward  and  backward  raster  scanning,  and  serpentine 
scanning. 

It  follows  from  the  definition  of  that  the  computation  of  a new  pixel  value  c(y)  can 
be  done  only  after  all  its  predecessors  (ordered  by  X)  have  been  computed.  Thus  recur- 
sive template  operations  may  not  be  performed  in  a globally  parallel  way  as  nonrecursive 
template  operations. 

Note  that  if  the  recursive  template  t is  defined  in  such  a way  that  5(t^y)  = 0,  the  empty 
set,  for  all  y E Y,  then  we  would  get  the  usual  nonrecursive  template  operations: 

a®-<t  = {(y,c(y)):c(y)=  a(x)  • t*y(x),  y € Y}. 

X€<5(t^y) 

So,  the  recursive  template  operations  are  natural  extensions  of  the  nonrecursive  ones. 

Similarly,  we  can  define  recursive  additive  maximum  and  recursive  multiplicative  maxi- 
mum as  follows.  If  a E R^  and  t E (R^,  rY^)^,  then 

a0^t  = {(y,  c(y)) : c(y)  = \f  a(x)  + t/y(x)  V \/  c(z)  + txy(z),  y G Y} 

xe5(t^y)  Z6«S(t^y) 

and  for  a £ (R^)x  and  t £ ((R^)x,  (R^)Y)Y, 

a©xt  = {(y,c(y)):c(y)=  \J  a(x)  • t/y(x)  V \J  c(z)  • txy(z),  y E Y}. 

x^5(t/y)  z65(t^y) 
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The  generalized  recursive  convolution  is  a linear  operation  while  the  other  two  are  non- 
linear operations. 

Note  that  we  can  generalize  the  above  three  basic  recursive  template  operations  by  defin- 
ing a generalized  recursive  product  between  images  and  recursive  templates  similar  to  the 
nonrecursive  case.  Let  7 be  an  associative  and  commutative  binary  operation  on  the  value 
set  F and  T the  global  reduce  operation  given  in  Chapter  2.  Let  Fi  and  F2  be  two  value  sets 
and  oj 9 o2  : Fi  x F2  — > F two  binary  operations.  If  a G Fx,  and  t G (FX,FY)Y,  then  the 
generalized  recursive  right  product  is  the  binary  operation  ©_<  : Fx  x (FX,FY)Y-+FY 
defined  by 

a©xt  = {(y,c(y)) : c(y)  = rx65(t7<y)a(x)  oa  t,<y(x)  7 rz65(t^y)c(z)  o2  txy(z),y  G Y}. 

Again,  in  addition  to  the  three  basic  recursive  template  operations,  a wide  variety  of 
recursive  template  operations  may  be  derived  from  this  generalized  recursive  template  prod- 
uct by  changing  the  binary  operations  7,0!,  and  o2.  We  will  present  some  applications  of 
the  generalized  recursive  template  operations  in  Chapter  6. 

3.3  Operations  between  Recursive  Templates 

Template  composition  and  decomposition  provide  a powerful  tool  in  the  image  algebra 
for  algorithm  optimization.  It  has  been  shown  [3,29]  that  for  a G RX  and  s,t  G (RX)Y, 

a©(s©t)  = (a©s)©t 

a©(s  + t)  = (a©s)  -f  (a©t) 
and  if  a G R^  and  s,t  G (RXOQ)Y, 

a0(sEt)  = (a0s)[V]t 


a[V](s  V t)  = (a[V]s)  v (a[V]t). 
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The  computation  of  (a©s)©t  is  usually  much  faster  than  the  direct  computation  of 
affi(s©t).  To  extend  the  template  decomposition  techniques  to  recursive  templates,  we 
need  to  define  the  operations  between  recursive  templates  and  prove  the  associative  and 
distributive  laws.  The  latter  will  be  done  in  Chapter  4. 

Let  X,  W C R71  be  finite.  Let  t £ (RX,RX)X,  s £ (R^,RX)X  be  two  recursive 
templates.  We  define  the  generalized  recursive  convolution  of  two  recursive  templates 
r = s©^t,  where  r = (r^,  r^)  £ (Rw,  Rx)x  , as  follows 


r-<  = 1 - (1  -8^0(1  -t«), 


where  1 £ (Rx)x  is  the  unit  template  and  is  defined  by 

Ix(x')  = 

We  also  define  the  addition  of  two  recursive  templates  r = s -f^/  t by 


1 if  x'  = x 
0 otherwise. 


r7<  — s/ ®(1  - *-<)  + (1  ~ s-<)®t7< 

= l-(l-s^)©(l-t^). 

Note  that  the  recursive  template  composition  is  defined  in  terms  of  nonrecursive  ones. 
Similarly,  if  t £ ((R|^)X,  (R^L)X)X,  s £ ((R$^)W,  (Rp^)X)X,  are  two  recursive 
templates  that  admit  the  partial  order  -<  on  the  point  set  X,  then  we  define  the  recursive 
template  composition  r = s©xt,  where  r = (r,<,r^)  £ ((R^)w,  (R^)X)X,  as  follows 

r*  = s^©t^ 

r ^ = O*A(lVsx)0(lVtx), 

where  0 £ is  the  zero  template  and  is  defined  by 

f 0 if  x'  = x 
1 — oo  otherwise. 


Ox(x')  = 
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Finally,  if  t G (R^,  R^)^,  s G ((R^,,  R^)*,  are  two  recursive  templates  that  ad- 
mit the  partial  order  -<  on  the  point  set  X,  then  we  define  the  recursive  template  composition 
r = s[V]-<t,  where  r = (r^,r^-)  G ((R^,,  R^)^-,  as  follows 

r*  = 

r~<  = <f>  A (0  V s_<)0(O  V t_<), 

where  (f>  G (R±oo)X  is  the  template  defined  by 

, / i\  _ \ if  x'  = x 

x X | +oo  otherwise. 

We  define  the  maximum  of  two  recursive  templates  r = s Vx  t by 


= 8^0(0  7^)7(0  7 8^)0^ 
= </>A(O7s^)0(O7t^). 


3.4  Computation  of  Recursive  Template  Operations 

The  computation  of  a recursive  operation  is  not  globally  parallel  because  the  order  of 
execution  is  restricted  by  the  partial  order  imposed  on  the  point  set  of  the  output  image. 
Note  that  any  poset  has  a linear  extension.  Thus,  a recursive  operation  may  be  computed 
in  any  one  of  these  linear  orders  sequentially  on  a conventional  computer.  In  the  case  that 
the  underlying  partial  order  defined  on  the  point  set  is  not  a total  or  linear  order,  some 
level  of  parallelism  can  be  achieved  by  using  parallel  architectures. 

Let  X be  an  n X n array.  Let  a G RX  and  ti,t2>l3  £ (R^)^  be  simplified  recursive 
templates,  whose  supports  are  given  below: 


tl 


w 
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Note  that  they  represent  three  typical  classes  of  recursive  templates,  namely,  one- 
dimensional, quarter-plane,  and  nonsymmetric  half-plane  recursive  templates.  Each  tem- 
plate t;  admits  a partial  order  (X,  which  can  be  defined  as  follows: 

1-  (X,-<tl):  (*',  (*,  j)  if  and  only  if  j'  < j , 

2-  (X,  -<t2):  (*/>J,)^t2(*"»j)  ^ and  only  if  one  of  the  following  conditions  holds 

(a)  i'  < i and  j'  < j , or 

(b)  i'  = i and  j'  < j , and 

3.  (X,-<t3):  {i',  j')~<t3(h  j)  if  and  only  if  one  of  the  following  conditions  holds 

(a)  i'  < i and  i'  + j'  < i + j,  or 

(b)  i'  = i and  j'  < j. 

Now  if  we  compute  bj  = a©X(.  ti  on  a sequential  machine,  then  it  may  be  realized 
on  any  linear  extension  of  (X,  ^tj)-  Of  course,  a natural  choice  would  be  forward  raster 
scanning.  The  total  time  required  to  compute  bj  will  be  0(n2). 

On  the  other  hand,  if  some  appropriate  parallel  architectures,  such  as  mesh  and  systolic 
arrays,  are  used,  then  the  total  running  time  would  be  linear,  which  can  be  shown  pictorially 
as  follows 


' 1 2 . 

1 2 . 

. n 

. n 

' 1 2 
2 3 

. n 

. n + 1 

' 1 2 . 

3 4 . 

. n 

. n + 2 

. 1 2 • 

. n 

n n + 1 . 

. 2n-l 

2n  — 1 2 n . 

. 3n  — 2 

where  the  value  in  each  entry  of  the  above  arrays  represents  the  time  that  the  new  pixel 
value  of  bj  at  that  location  may  be  computed. 

A new  pixel  can  be  processed  once  all  the  pixels  ordered  previous  to  it  by  the  given 
partial  order  have  been  processed.  Those  pixels  that  are  not  comparable  under  the  given 
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partial  order  do  not  depend  on  each  other;  therefore,  they  can  be  processed  simultaneously 
on  parallel  architecture.  Thus,  it  requires  n unit  times  to  compute  bi,  2n  - 1 unit  times 
to  compute  b2,  and  3n  — 2 unit  times  to  compute  b3.  Here,  by  one  unit  time  we  mean  the 
time  to  compute  a single  new  pixel  value  of  the  output  image. 


CHAPTER  4 

ALGEBRAIC  PROPERTIES  OF  RECURSIVE  OPERATIONS 


In  Chapter  3,  we  have  defined  the  recursive  templates  and  recursive  template  operations. 
In  this  chapter,  we  will  investigate  their  mathematical  structures  and  algebraic  properties. 
Some  of  these  algebraic  properties  are  essential  to  develop  template  decomposition  tech- 
niques for  algorithm  optimization  and  parallel  image  processing. 

4.1  Linear  Recursive  Operations 

Linear  transforms  are  the  most  understood  and  commonly  used  operations  in  image 
processing.  In  this  section,  we  will  show  the  relationship  between  the  linear  recursive 
operations  in  the  image  algebra  and  matrix  operations  in  linear  algebra. 

It  is  known  that  the  image  algebra  has  a close  relationship  with  linear  algebra.  In 
particular,  Gader  [28]  showed  that  the  structure  ((F^)^,  ©,  +)  is  a ring  and  it  is  isomorphic 
to  the  ring  structure  (Mn,  *,  +)  in  linear  algebra,  where  Mn  is  the  set  of  all  n X n matrices 
with  values  in  F,  and  n = card(X),  the  cardinality  of  X.  Since  the  recursive  template  and 
generalized  recursive  convolution  are  more  general  than  their  non-recursive  counterparts, 
one  would  expect  that  only  some  subclasses  of  recursive  templates  could  enjoy  such  an 
algebraic  structure.  In  this  section,  we  will  show  that  this  is  indeed  the  case.  The  set  of 
invariant  quarter-plane  recursive  templates  equipped  with  recursive  template  convolution 
and  addition  forms  a commutative  ring.  This  special  class  of  recursive  templates  corresponds 
to  the  class  of  quarter-plane  recursive  filters. 

Let  X be  a partially  ordered  point  set  with  the  partial  order  -<,  and  F a value  set.  Let 
n — card(X).  Note  that  if  -4  is  an  empty  partial  order,  i.e.,  no  element  is  comparable  to 
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any  other  element  in  X,  then  X would  simply  be  a regular  point  set.  We  know  that  every 
partial  order  has  a linear  extension.  Let  xa,  x2, . . . , xn  be  a linear  extension  of  -<  on  X.  We 
define  a mapping  v : Fx  — *■  F by  i/(a)  = (a(xi),  a(x2), . . . , a(x„)).  For  our  convenience, 
we  will  consider  v{&)  as  a row  vector  instead  of  column  vector.  Clearly,  v is  1-1  and  onto.  It 
maps  a set  of  images  Fx  to  a set  of  n— dimensional  vectors.  Note  that  given  a partial  order 
-<,  there  can  be  many  different  linear  extensions  of  -<.  Hence,  the  corresponding  mapping 
v can  be  different  from  one  to  another  by  a permutation. 

For  each  template  t 6 (Fx)x,  we  define  a matrix  Mt  = (m,j),  where  m,y  = tXj.(xj). 
Note  that  the  jth  column  of  Mt  is  (^(tx^))'.  Let  Mn  be  the  set  of  n x n matrices  with  values 
in  F.  Define  a mapping  $ : (Fx)x  -►  Mn  by  $(t)  = Mt.  Gader  has  shown  the  following 
theorem  [28]: 

Theorem  4.1  (Gader)  is  a ring  isomorphism  o/((Fx)x,  ©,+)  onto  (Mn,*,+). 

We  will  show  a similar  result  for  recursive  templates  and  operations. 

Definition  An  n x n matrix  A = (d,j)  is  said  to  be  upper  triangular  if  a,j  = 0 for  all 
i > j,  strictly  upper  triangular  if  atJ  = 0 for  all  i > j,  and  unit  upper  triangular  if  at]  = 0 
for  all  i > j and  au-  = 1 for  all  i. 

Definition  An  n x n matrix  A = (a,y)  is  said  to  be  a Toeplitz  matrix  if 


«0 

G1 

an-2 

®n— 1 

a-i 

do 

®n— 2 

a-n+2 

. . . 

do 

dl 

. n+1 

a-n+2  ■ • • 

d_l 

Oq 

That  is,  aij  = a,-+fc,j+jt  for  all  0 < i,j,k  < n - 1 where  Q<i-\-k<n-l  and 
0<j  + k<n  — 1.  A Toeplitz  matrix  is  constant  along  the  diagonals. 
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Definition  An  n x n matrix  A = (a,j)  is  said  to  be  a block  Toeplitz  matrix  with  Toeplitz 
blocks  if 

Ao  Ax  • • • Am—2  -“dm— 1 

A-x  Aq  ...  Am—2 


A = 


5 


• • • Ao  A\ 

. . . A_x  A0 

where  A*  is  a Toeplitz  matrix  for  all  — m + 1 < k < m — 1. 

A is  said  to  be  an  upper  triangular  block  Toeplitz  matrix  with  upper  triangular  Toeplitz 

blocks  if 


A-m+2 

A— m+l  A_ m+2 


Ao 


Ax 

Aq 


Am— 2 Am_x 
Am— 2 


Ao  Ax 
Ao 

where  A*  is  an  upper  triangular  Toeplitz  matrix  for  all  0 < k < m - 1. 

Let  be  the  set  of  all  n X n upper  triangular  Toeplitz  matrices.  Let  Tn  be  the  set  of 
all  n x n upper  triangular  block  Toeplitz  matrices  with  upper  triangular  Toeplitz  blocks, 
where  each  block  is  of  the  size  \/n  X y/n.  We  assume  that  these  matrices  are  defined  on  the 
value  set  R. 

The  following  theorem  is  given  in  Gohberg  et  al.  [39]. 


Theorem  4.2  The  set  T'n  is  an  algebra,  that  is,  it  is  closed  under  the  operations  of  addition, 
multiplication  by  scalars,  and  matrix  multiplication.  Moreover,  if  A £ T'n  and  det(A)  / 0, 
then  A-1  € T'. 

Corollary  4.3  is  a commutative  ring. 

Proof:  Let  A,  B £ T„,  C = AB,  and  D = BA. 


ao  a\  . 

• 2 

Q"n—  1 

bo  bi  . 

• bn-2 

bn- 1 

a0  . 

• 

^n—2 

bo  . 

bn- 2 

ao 

ax 

a0 

B = 

bo 

h 

bo 
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By  Theorem  4.2,  both  C and  D are  in  T'n.  Thus,  we  need  only  to  check  if  coj  = doj  for 
0 < j < n — 1.  Note  that 

n— 1 

COj  = ^2  aOkbkj 

k= 0 

= aobj  + aibj-i  + b djbo 

n— 1 

doj  — ^ ] bokakj 
k= o 

= b0aj  + &1U/-1  + • • • + bjdQ. 

Hence,  we  have  c0j  = d0j  for  0 < j < n — 1. 

Q.E.D. 

From  this,  we  can  derive  the  following  theorem,  which  we  are  going  to  use  later  in  this 
section. 

Theorem  4.4  (Tn,*,+)  is  a commutative  ring. 

Proof:  Let  A,BeTn  and  C = AB.  Then,  for  * < j , 

n— 1 

Ctj  = AikBkj 

k= 0 

= AoBj-i  + + • • • + Aj^iBg. 

Since  A,  B £ Tn,  C,+/J+/  = C4j  and  C,y  = 0 for  i > j.  Since  Aik,Bkj  6 T£,  by  Corollary 
4.3,  C,j  is  in  T^.  Thus,  C is  in  Tn,  and  so  Tn  is  closed  under  the  matrix  multiplication. 
From  the  above  equation  and  the  fact  that  the  matrix  multiplication  is  commutative  in  T'n, 
we  know  that  the  matrix  multiplication  is  commutative  in  Tn  also. 

It  is  easy  to  show  that  Tn  is  closed  under  matrix  addition  and  it  forms  a ring  with  matrix 
addition  and  multiplication. 


Q.E.D. 
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Definition  An  invariant  template  t £ is  said  to  be  quarter-plane  if  its  support 

S(to)  lies  in  a quarter  plane,  where  0 = (0,0).  In  particular,  t is  said  to  be  causal  if  its 
support  lies  in  the  first  quarter  plane.  We  denote  the  set  of  all  causal  templates  in  (Rx)x 

by  TX- 

The  following  theorem  is  from  Gader  [28], 

Theorem  4.5  If  t G (Fx)x  be  translation  invariant  and  Mt  be  the  matrix  corresponding 
to  t,  that  is,  Mt  = 'P(t),  then  Mt  is  block  Toeplitz  with  Toeplitz  blocks. 

Corollary  4.6  If  t £ Tx  and  Mt  = (m,y)  be  the  matrix  corresponding  to  t,  then  Mt  is 
upper  triangular  block  Toeplitz  with  upper  triangular  Toeplitz  blocks,  assuming  that  we  order 
X in  the  forward  raster  scanning  order. 

Proof:  By  theorem  4.5,  Mt  is  block  Toeplitz  with  Toeplitz  blocks.  Let  Mt  = (Aty),  0 < 
i,j  < n - 1,  where  each  block  A.y  = ( aki ) is  an  y/n  X y/n  matrix.  Note  that  the  rth  column 
of  Mt  is  (^(tXr))'.  If  we  order  X in  the  row  major  form,  then  a^  = TOtn+fcjn+;  = t^(i,k). 
The  causality  of  t implies  that  t (j,i)(i,k)  ^ 0 if  and  only  if  i < j and  k < l.  It  follows  that 
A{j  = 0 for  i > j and  a m — 0 for  k > / within  each  block  A,-j . Thus,  Mt  is  upper  triangular 
block  Toeplitz  with  upper  triangular  Toeplitz  blocks. 

Q.E.D. 

The  following  is  immediate  from  the  above  theorems. 

Theorem  4.7  \I>  is  a ring  isomorphism  o/(Tx,  ©,+)  onto  (T„,*,+). 

We  will  use  these  theorems  to  derive  the  structure  for  recursive  templates  and  operations. 
Let  Pn  be  the  set  of  ordered  pairs  (A,  B)  where  A,B  £ Tn  such  that  det(B)  ^ 0.  We 
define  a relation  ~ on  Pn  as  follows: 


(A,fl)~(C,D) 


if  and  only  if  AD  = BC 
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We  claim  that  this  defines  an  equivalence  relation  on  Pn. 

1.  If  (A,  B ) 6 Pn,  then  (A,  B ) ~ (A,  B)  since  AB  = BA  by  Corollary  4.4. 

2.  If  (A,  B),(C,  D)  £ Pn  and  ( A,B ) ~ (C,  D),  then  AD  = BC.  Since  the  matrix 
multiplication  in  Tn  is  commutative,  we  have  CB  = DA,  and  so  ( C , D)  I *'w»  (■ A,B ). 

3.  If  (A,  B ),  {C,  D),  (E,  F ) G Pn  and  (A,  B)  ~ (C,  D ) and  (C,  D)  ~ (E,  F ),  then  AD  = 
BC  and  CF  = DE,  hence  BCF  = BDE.  From  AD  = BC  and  BCF  = BDE,  we 
have  ADF  = BDE.  Since  the  matrix  multiplication  in  Tn  is  commutative  and  D is 
invertible,  it  follows  that  AF  = BE.  Thus,  (A,  B ) ( E,F ). 

By  showing  the  reflexivity,  symmetry,  and  transitivity  of  the  relation  we  have  proven 
that  it  is  indeed  an  equivalence  relation  on  Pn. 

Let  [A,  B]  be  the  equivalence  class  of  (A,  B)  in  Pn,  and  let  Pn  be  the  set  of  all  such 
equivalence  classes  [ A,  B ] where  ( A,  B ) G Pn ■ We  define  an  addition  and  a multiplication 
over  the  elements  in  Pn  as  follows: 

[A,  B]  + [C,  D]  = [AD  + BC,  BD] 

[A,  B]  * [C,  D]  = [AC,  BD] 

Since  det(B)  / 0 and  det(D)  / 0,  det(BD)  ± 0,  and  so  both  [AC,  BD]  and  [AD  + BC,BD] 
are  in  Pn . 

Next,  we  show  that  the  addition  is  well  defined,  that  is,  if  [A,  B]  = [A',  B']  and  [C,  D]  = 
[C' , D'],  then 

[A,B]  + [C,D]  = [A',B']A[C',D']. 

Since  [A,  B]  = [A',B']  and  [C,  D]  = [C',D'],  we  have  that  AB1  = BA’  and  CD'  = DC'. 
It  follows  that  ADB'D'  = BDA'D'  and  BCB'D'  = BDB'C'.  Thus,  (AD  + BC)B'D'  = 
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BD(A'D'  + B'C' ),  which  implies  that  [AD  + BC,BD]  = [A'D1  + B'C',B'D'],  that  is, 
[A,B]  + [C,D]  = [Af,B']  + [Ct,D’\. 

Similarly,  we  can  show  that  the  multiplication  is  well  defined,  that  is,  if  [A,  B]  — [A',  B '] 
and  [C,  D]  = [C , D1],  then 

[A,B]*[C,D]  = [A\B'}*[C',D'). 

Since  [A,  B]  — [A',B']  and  [C,  D]  = [C',D'],  we  have  that  AB'  = BA'  and  CD 1 = DC1. 
Thus,  ACB'D'  = BDA'C',  which  implies  that  [AC,  BD]  = [A'C' ,B' D'],  that  is,  [A,  B]  * 
[C,D]  = [A',B']*[C',D'). 

The  following  theorem  says  that  the  set  Pn  together  with  the  addition  and  multiplication 
defined  above  forms  a commutative  ring. 

Theorem  4.8  (P„,*,+)  is  a commutative  ring  with  identity. 

Proof: 

Let  [A,B],[C,D],[E,F]e  Pn.  First,  we  show  that  ( Pn , +)  is  an  abelian  group. 

1.  Since  Tn  is  closed  under  matrix  addition  and  multiplication  by  Theorem  4.4,  both 
AD  + BC  and  BD  are  in  Tn.  Furthermore,  det(B)  ± 0 and  det(D)  / 0 imply  that 
det(BD)  ± 0.  Thus  [A,  B]  + [C,  D]  = [AD  + BC , BD]  is  in  Pn. 

2.  Associativity 

([A,B]  + [C,D])  + [E,F]  = [ AD  + BC,BD]  + [E,F ] 

= [(AD  + BC)F  + BDE,  BDF] 

= [ADF  + BCF  + BDE,  BDF] 

[A,5]  + ([C,D]  + [£,F])  = [A,B]  + [CF  + DE,DF] 

= [ADF  + B(CF  + DE),  BDF] 
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= [ADF  + BCF  + BDE,  BDF] 

= ([  A,B]  + [C,D])  + [E,F\. 

3.  By  Theorem  4.4,  Tn  is  a commutative  ring.  So  we  have 

[A,B\  + [C,D\  = [AD  + BC,  BD] 

= [CB  + DA,DB] 

= [C,D]+[A,B]. 

4.  Let  [0,  G ] G Pn.  Since  (0,  G)  ~ (0,  G ')  for  any  G and  G',  [0,  G]  is  unique  in  Pn.  We 
denote  it  by  [0,/n],  where  In  is  the  identity  matrix  in  Tn. 

[A,B]  + [0,In]  = [ AIn,BIn ] 

= [A,B\ 

[0  ,In\  + [A,B]  = [InA,InB] 

= [A,B]. 

Thus,  [0,/n]  is  the  identity  for  the  operation  +. 

5.  Inverse 


[A,B]  + [-A,B]  = 


[ -A,B]  + [A,B ] = 


[AB  + B(-A),  BB] 
[0  ,BB] 

[0  ,/„] 

[-AB  + BA,BB] 

[0  ,BB] 

[0,4]. 


Thus,  the  additive  inverse  of  any  [A,B]  in  Pn  exists  and  it  is  [-A,  B], 
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Next,  we  need  to  show  that  ( Pn , *)  is  a semi-group  with  identity. 

Clearly,  Pn  is  closed  under  *.  Note  that 

([A,B]*[C,D])*[E,F\  = [ AC,BD]*[E,F ] 

= [ACE,BDF] 

= [A,  B]  * [ CE , DF] 

= [A,B]*([C,D}*[E,F]). 

Thus,  the  operation  * is  associative.  Note  that  [/„,/„]  is  in  Pn  and 

[A,B]*[In,In\  = [ A,B } 

= [In->In\*\A,B\. 

Therefore,  [/„,/„]  is  the  identity  for  the  operation  *. 

Finally,  we  show  that  the  operation  * is  distributive  over  the  operation  +. 

[A,B]*([C,D]  + [E,F])  = [A,B]*[CF  + DE,DF] 

= [ACF  + ADE,  BDF] 

[A,B]*[C,D]  + [A,B]*[E,F]  = [AC,  BD)  + [AE,  BF] 

= [ACBF  + BDAE,  BDBF] 

= [ACF  + DAE,  BDF } 

= [A,B}*([C,D]  + [E,F]). 

Similarly,  it  is  easy  to  verify  the  right  distributive  law. 

Q.E.D. 

Corollary  4.9  ( Pn , *,  +)  is  a field,  where  Pn  is  a subset  of  Pn  such  that  for  any  [A,  B ] E Pn, 
det(B)  0,  and  det(A ) — 0 if  and  only  if  A — 0. 
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Proof: 

If  [A,  B]  £ Pn  and  [A,  B]  ± [0,  Jn],  then  det(A)  0,  and  so  [B,A]  £ Pn.  Also, 

[A,  B]  * [U,  A]  = [AB,  BA] 

= [B,A]  * [A,B]. 

Thus,  the  multiplicative  inverse  of  any  non-zero  element  [A,  B]  exists  and  it  is  [B,A\. 

Q.E.D. 

Definition  An  invariant  recursive  template  t £ (FX,FY)^  is  said  to  be  quarter-plane  if 
both  <S(t^o)  and  <S(txo)  he  in  the  same  quarter-plane,  where  0 = (0,0).  In  particular,  t 
is  said  to  be  causal  if  both  and  are  causal.  We  denote  the  set  of  all  causal  recursive 
templates  by  T^,x. 

Define  a mapping  $ : TxX  Pn  by  $(t)  = ($(t ^),In  - $(t^)).  Also,  define  a relation 
- on  TxX  as  follows: 

s ~ t if  and  only  if  3>(s)  ~ $(t) 

Since  ~ is  an  equivalence  relation,  the  relation  ~ is  an  equivalence  relation  on  T^x.  We 
denote  the  set  of  all  equivalence  classes  over  TxX  by  TxX. 

Define  a mapping  $ : TxX  -*■  Pn  by  $([tj)  = [*(t/),/n-®(tx)].  The  following  theorem 
characterizes  the  algebraic  structure  of  a class  of  causal  recursive  templates. 

Theorem  4.10  4>  is  a ring  isomorphism  of  (T^x,  onto  (Pn,*,+). 

Proof: 

If  [s],[t]  £ T^x,  then 

+-<W)  = *([H+«t]) 
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= *([(«*©(  1 - t«)  + (1  - s^)©^,  1 - (1  - sx)©(l  - t^))]) 

= [*(8/®(i  - u)  + (i  - s,,)®^), /„  - 9(1  - (i  - s^)®(i  - u))]. 
By  Theorem  4.7,  we  have 

Hi*}  W)  = m*A)Vn  ~ *(t^))  + (In  ~ ¥(SX))* (trf),  (In  ~ 9(aJ)(In  - 9(t<))] 

= [*M,(in  - ns*))}  + muWn  - nu))] 

= $([s])+$([t]). 


Similarly,  we  have 


*([s]©-c[t])  = 4>([s®xt]) 

= ^([(S,<©t/,  1 - (1  - sx)©(l  - t^))]) 

= [$(s/®t^),/n  - 9(1  - (1  - SX)©(1  - tx))] 
= [9(at)9(tt),  (In  - 9(a  J)(In  - *(t^))] 

= [tf  (s,<),  (In  - W(sx))]  * [*(t^),  (In  - *(tX))] 

= *([s])*$([t]). 


Let  [A,  B]  G Pn • Since  [AX,  BX]  = [A,B\  and  det(B)  0,  without  loss  of  generality, 
we  may  assume  that  the  elements  of  the  main  diagonal  of  B are  all  Ts.  Rewrite  [A,  B] 
as  [A,  In  - C],  then  C is  strictly  upper  triangular.  It  follows  that  t = (\E,_1(A),  $_1(C)) 
defines  a causal  recursive  template  in  TxX,  where  4>(t)  = [A,/n  - C],  and  so  $ is  onto. 

Suppose  that  $([s])  = $([t]).  Then,  $(s)  ~ $(t).  By  definition,  we  will  have  s ~ t, 
which  implies  that  [s]  = [t].  Thus,  4>  is  1-1  and  onto. 

The  unit  recursive  template  in  TxX,  denoted  by  1^,  is  given  by  1^  = (1,0),  where  0 
is  the  null  or  empty  template  in  (Rx)x  (the  support  of  0 is  an  empty  set).  By  definition, 
$([!-<])  = [In,In],  which  is  the  unit  element  in  Pn. 


Q.E.D. 
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Corollary  4.11  4>  is  a field  isomorphism  of  { TxX,  ©-<,+-<)  onto  (Pn,*,+),  where  T^x 
is  a subset  of  TxX  such  that  for  any  t £ T^x,  except  0X  = (0,0),  t^x(x)  / 0,  Vx  £ X. 

Now  the  mathematical  structure  of  the  linear  recursive  templates  and  operations  is  clear. 
We  can  use  this  result  to  show  some  useful  properties  of  recursive  template  operations.  Since 
(T^x,  ©^,  is  a commutative  ring,  we  have 

r©_<(s©_<t)  = (rffixs)©^t 

s©^t  = t ©^/S 

r©_<(s  t)  = (r©^s)  (rffi^t) 

for  any  r,s,t  £ T^X. 

Theorem  4.12  If  a £ Rx  and  t £ TxX,  then 

Kaffi^t)  = i/(a)«(t/)(In  - 

Proof: 

If  b = a©^t,  then 

z/(b)  = K*W/)  + Kb)^(t^). 

Since  t £ TxX,  $(t^)  is  strictly  upper  triangular,  having  all  zeros  along  its  main  diagonal, 
and  so  In  — KK)  is  invertible.  Thus,  we  have 

Kb)  = Ka)*(t /)(/„  - ®(t^))-1. 


Theorem  4.13  If  a £ Rx  and  s,  t 6 TxX,  then 


Q.E.D. 


a©^(sffi^t)  = (affi^s)©^t. 
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Proof: 

Note  that  s©^t  = (s^  ©t^,  1 — (1  — s^)©(l  — t^)).  By  Theorem  4.12  and  Theorem  4.7, 

Ka®-<(s®-<t))  = Ka)«r(s7<©t7<)[/n-$(l-(l-Sx)ffi(l-tx))]-1 
= i/(a)*(s/)*(t/)[(/B  - *(s^))(/n  - 
= u(a )9M(In  - *(s^))_1*(t^)(/n  - 
= Ka®^s)^(t^)(/n  - ^(t^))-1 
= i/((affi^s)©^t). 


Q.E.D. 

Theorem  4.14  If  a £ Rx  and  s,  t e T^x,  then 

a©^(s  +x  t)  = (affixs)  + (a©^t). 

Proof: 

First,  note  that  s t = (s^©(l  - t^)  + (1  - s^)©t^,  1 — (1  - s^)©(l  - t^)).  By  The- 
orem 4.12  and  Theorem  4.7, 

i/(a©^(s  +,.  t)) 

= ^^©(l  - U)  + (1  - 8^)©t/)[/n  - *(l  - (1  - S^)ffi(l  - t^))]-1 

= Ka)(*M(/n  ~ *(t„))  + (In  - *(B^))*(t/))[(/n  - *(S ■<))(/„  - * (t^))]-1 
= */(a)«(B *)(Jn  - tf(s^))-1  + i/(a)¥(t,0(JB  - tf(t^))-1 
= i^(a©^s  -f-  a©xt). 


Theorem  4.15  If  a,  b G RX  and  t G T^,x>  then 


Q.E.D. 


(a  + b)©^-t)  = (a©^t)  + (b©^t) 
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Proof: 

By  Theorem  4.12, 

i/((a  + b)©xt)  = i/(a  + b)$(t/)[JB-*(t^)]-1 

= (i/(a)  + i/(b))*(t/)[/B-*(t^)]-1 
= ^(affi^t  + b©^t). 

Q.E.D. 

By  Corollary  4.11,  each  causal  recursive  template  t has  an  inverse  t-1  if  t is  not  an 
empty  template  and  t^x(x)  ^ 0.  Note  that 

Since  t^x(x)  ^ 0,  the  inverse  of  4>([t])  exists.  Multiply  ($([t]))_1  on  each  side  of  the 
equation  above,  we  have 

^([t-1])  = ($([t]))_1  * ^([tP^t-1]) 

= [/n-¥(tx),*(t,<)]. 

Thus, 

[t-1]  = *-\[in-*(u)MU)}) 

= [(l-t^l-trf)]. 

Note  that  a causal  nonrecursive  template  t with  tx(x)  ^ 0 is  always  invertible  since 
^(t)  is  invertible.  It  can  be  viewed  as  a causal  recursive  template  (t^,0).  If  t is  invertible, 
then  there  exist  a causal  nonrecursive  template  ^ and  a causal  recursive  template  t2  such 
that  if  we  consider  t as  a nonrecursive  template,  then 


t©ti  = 1 


43 


and  if  we  consider  t as  a recursive  template,  then 

t®«t2  = 1^ 

Therefore,  we  have  two  ways  to  compute  the  inverse  of  a causal  nonrecursive  template 


1.  t"1  = ^((^(t))-1) 

2.  t-1  = (1-^,1-^) 

such  that  either  t®t-1  = 1 or  t®xt-1  = lx. 

The  first  method  gives  a causal  nonrecursive  template  and  the  second  one  yields  a causal 
recursive  template.  In  general,  the  size  of  the  support  of  the  inverse  template  generated  by 
the  second  method  will  be  the  same  as  the  given  nonrecursive  template,  whereas  the  the 
size  of  the  support  of  the  inverse  template  generated  by  the  first  method  will  be  as  large  as 
the  underlying  point  set  X. 

Example  Let  X be  a 3 x 3 point  set.  Let  t E (Rx)x  be  a causal  nonrecursive  template 
defined  as  follows: 
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and 


W1 


1 
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So,  the  inverse  template  would  be 


.375 

.375 
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.5 

.5 
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.5 

w 

Applying  the  second  method,  we  will  have  a simplified  recursive  inverse  template  as 


follows 


5 

5 

w 

Clearly,  the  advantage  of  considering  a causal  nonrecursive  template  t as  a causal  re- 
cursive template  and  computing  its  inverse  as  a causal  recursive  template  is  that  we  can 
get  a small  inverse  template  so  that  the  convolution  with  the  inverse  template  can  be  very 
efficient.  This  idea  has  been  widely  used  in  image  coding  and  decoding  techniques  [32,40], 
Before  we  leave  this  section,  we  would  like  to  show  that  although  the  theorems  we  have 
proven  in  this  section  are  related  to  causal  recursive  templates,  they  are  also  valid  for  any 
class  of  quarter-plane  templates. 
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Let  X C Z2  be  a 3 x 3 point  set,  and  s,t  G (Rx,  R^)^  be  simplified  recursive  templates 
defined  as  follows 
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3 
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4 

2 

The  matrix  representation  of  these  templates  would  be 
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The  transposes  of  s and  t and  their  matrix  representations  are 
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where  ^(s^)  and  ^(t^)  are  the  identity  matrix  /9. 
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Note  that  all  of  these  matrices  are  upper  or  lower  triangular  block  Toeplitz  matrices.  In 
general,  the  matrices  and  ^(t^)  corresponding  to  a recursive  quarter-plane  template 

t € (FX,FX)X  fall  into  one  of  the  following  cases: 

1.  upper  triangular  block  Toeplitz  with  upper  triangular  Toeplitz  blocks  if  the  supports 
of  t^o  and  t_<o  lie  in  the  first  quarter  plane 

2.  lower  triangular  block  Toeplitz  with  upper  triangular  Toeplitz  blocks  if  the  supports 
of  t^o  and  t_/o  lie  in  the  second  quarter  plane 

3.  lower  triangular  block  Toeplitz  with  lower  triangular  Toeplitz  blocks  if  the  supports 
of  t^o  and  t^o  lie  in  the  third  quarter  plane 

4.  upper  triangular  block  Toeplitz  with  lower  triangular  Toeplitz  blocks  if  the  supports 
of  t^o  and  txo  lie  in  the  fourth  quarter  plane 

It  is  easy  to  see  that  if  we  let  Tn  be  any  one  of  these  classes  of  matrices,  all  the  theorems 
would  be  true  as  well.  That  is,  the  results  in  this  section  apply  to  any  class  of  quarter-plane 
templates. 

4.2  Non-linear  Recursive  Operations 

In  this  section,  we  will  show  the  relationship  between  the  non-linear  recursive  operations 
in  the  image  algebra  and  the  matrix  operations  in  the  minimax  algebra. 

The  relationship  between  the  image  algebra  and  the  minimax  algebra  was  established 
by  Davidson  in  her  dissertation  [29].  In  particular,  let  (F,  x,  V)  and  (F,  *,  V)  be  two  belts. 
Let  $ be  the  same  mapping  defined  in  Section  4.1.  Davidson  proved  the  following  theorems: 

Theorem  4.16  $ is  a belt  isomorphism  of  ( Fx)x,  0,  V)  onto  (M„,+,  V). 

Theorem  4.17  $ is  a belt  isomorphism  of  ( Fx)x,  ©,  V)  onto  (Mn,*,v). 
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where  Mn  is  the  set  of  all  n X n matrices  with  values  in  the  set  F,  and  n = card(X). 

We  will  extend  this  result  to  non-linear  recursive  operations.  We  show  that  the  set 
of  invariant  quarter-plane  recursive  templates  equipped  with  recursive  template  operations 
additive  maximum  and  maximum  forms  a commutative  belt. 

Let  (R_oo,  X,  V)  be  a belt  with  the  zero  element  — oo  and  the  unit  element  0.  So,  the 
operation  X is  really  the  regular  addition  +.  Let  T'n  be  the  set  of  all  n X n upper  triangular 
Toeplitz  matrices  with  values  in  the  set  R-qq.  Let  Tn  be  the  set  of  all  n x n upper  triangular 
block  Toeplitz  matrices  with  upper  triangular  Toeplitz  blocks  with  values  in  the  set  R_0O, 
where  each  block  is  of  the  size  y/n  X y/n. 

Theorem  4.18  (T',  X,V)  is  a commutative  belt. 

Corollary  4.19  (Tn,  x,V)  is  a commutative  belt. 

Let  0n  and  On  denote  the  zero  and  unit  elements  in  (Tn,  X,  V),  respectively.  Then,  0n 
is  the  matrix  with  all  entries  being  — oo  and  the  unit  element  On  is  the  matrix  with  0’s  on 
its  main  diagonal  and  — oo  elsewhere.  Let  Tx  be  the  set  of  causal  templates  in  (R^f- 

Theorem  4.20  'k  is  a belt  isomorphism  of  ( Tx,  El,  V)  onto  (Tn,  X,  V). 

The  proofs  of  the  three  theorems  above  are  similar  to  the  proofs  of  Theorem  4.3,  Theorem 
4.4,  and  Theorem  4.5,  and  thus  omitted  here. 

We  will  use  these  theorems  to  derive  the  structure  for  non-linear  recursive  templates 
and  operations. 

Definition  Let  A 6 Mn.  The  transitive  closure  of  A,  denoted  by  A+,  is  defined  as 

A+  - On  V A V A2  V • • • 

where  Ak  = A X A x ■ ■ ■ X A.  A is  said  to  be  stable  if  Ak  = 0n  for  some  finite  k. 

v ' ' 

k 

Note  that  if  A is  stable,  then  its  transitive  closure  exists  and  can  be  computed  in  finite 
number  of  matrix  multiplications  and  additions. 
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Lemma  4.21  If  A 6 Tn  and  the  main  diagonal  elements  of  A be  — oo,  then 

A+  = On  V A V A2  V • • • V A”-1 . 

Proof: 

Since  A E Tn,  A is  upper  triangular.  With  the  fact  that  the  main  diagonal  elements  of  A 
are  -oo,  it  follows  that  An  = 0n,  and  so 

A+  = On  V A V A2  V • • • V An-1 . 


Q.E.D. 

Let  Pn  be  the  set  of  ordered  pairs  ( A,B ),  where  A,  B E Tn  and  the  main  diagonal 
elements  of  B are  — oo.  We  define  a relation  ~ on  Pn  as  follows: 


(A,  B ) ~ (C,  D)  if  and  only  if  A x B+  = Cx  D+ 


Clearly,  this  is  an  equivalence  relation  on  Pn. 

Let  [A,  B]  be  the  equivalence  class  of  (A,B)  in  Pn,  and  let  Pn  be  the  set  of  all  such 
equivalence  classes  [A,  B]  where  ( A,  B ) E Pn.  We  define  two  operations  on  the  elements  in 
Pn  as  follows: 


[A,  B]  V [C,  D]  = [(A  X (On  V D ))  V ((On  V B)  X C),  <j>n  A ((0„  V B)  X ( On  V D))] 
[A,  B]  x [C,  fi]  = [Ax  C,  0n  A ((On  V P)  x (On  V £>))], 


where  (j)n  is  the  matrix  with  — oo  on  its  main  diagonal  and  +oo  elsewhere,  i.e., 


— oo  +oo  . . . +oo 
+oo  — oo  . . . +oo 

0 n — , . . 

+oo  +oo  ...  — oo 

Since  the  diagonal  elements  of  B and  D are  — oo,  we  can  see  that  so  is  the  matrix 


< K A ((0„  Vfl)x  (On  V £>)).  Thus,  the  set  Pn  is  closed  under  the  operations  V and  X defined 
above.  Also,  it  is  easy  to  show  that  the  operations  V and  X are  well  defined  on  Pn. 
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The  following  theorem  says  that  the  set  Pn  together  with  the  operations  V and  X defined 
above  forms  a commutative  belt. 

Theorem  4.22  (Pn,x,V)  is  a commutative  belt  with  identity. 

Proof: 

Let  [A,  F],  [C,  F],  [ E , F]  € Pn.  First,  we  show  that  ( Pn , V)  is  a commutative  band. 

1.  Associativity 

Note  that  On  V (<f>n  A ((On  V B)  x (Ob  V D)))  = (0B  V B)  x (On  V F).  Thus, 
([A,F]v[C,F])v[F,F] 

= [Ax  (0B  V F)  V (On  V B)xC,<t>n  A ((On  V 5)  X (On  V Z?))]  V [£,  F] 

= [(A  X ( On  vfi)v  (0„  V f?)  X C)  X (On  V F)  V (On  V <t>n  A ((On  V F)  X 
(On  V F)))  x£,4  A (On  V <t>n  A ((On  V F)  X (On  V F)))  X (On  V F)] 

= [Ax  ( On  V F)  X (On  V F)  V (On  V F)  X C X (On  V F)  V (On  v F)  X 
(On  V D)  X E,<t>n  A ((On  V F)  X (On  V F)  X (On  V F))] 

and 

[A,F]V([C,F]V[F,F]) 

= [A,  F]  V [C  x (On  V F)  v (On  V F)  x F,  <fin  A ((On  v F)  x (0B  v F))] 

= [A  X (On  v <j>n  A ((On  v F)  X (On  v F)))  V (0B  V F)  X (C  X (On  V F)  V 
(On  V F)  X F),  <t>n  A ((On  v F)  X (On  V 4>n  A ((On  V F)  X (On  V F))))] 

= [Ax  (On  V F)  X (On  VF)V(OnVF)xCX  (On  v F)  V (On  v F)  X 
(On  V F)  X F,  (t>n  A ((On  V F)  X (On  V F)  X (On  v F))] 


([A,F]V[C,F])V[F,F], 
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2.  Commutativity 

By  Theorem  4.19,  Tn  is  a commutative  belt.  So  we  have 


[A  X (On  V O)  V (On  V B)  X C,  <fin  A ((On  V B)  X (On  V 2?))] 
[C  X (On  Vfl)v  (On  V D)  X A,  <f>n  A ((0„  vD)x  (On  V £))] 
[C,0]v[A,0]. 

3.  Idempotency 
Note  that 


[A,  5]  V [0,0]  = 


A X (0„  V £)  X [<t>n  A (On  V 5)  X (0„  V £)]+ 
= A X (On  V J9)  X [<fin  A (On  V £ V B X £)]  + 

= Ax(OnV  B)x(Bv  B x B)+ 

= Ax(OnV  B)x  B+ 

= A x B+ . 


So,  by  the  definition  of  the  equivalence  relation  we  have 


[A  x (0„  V B ),  <fin  a ((On  V B)  x (0„  V 5))]  = [A,  B]. 

Thus, 

[A,0]V[A,0] 

= [A  X (On  V5)v  (On  V B)  X A,  (fin  A ((0„  VB)x  (0„  V 5))] 
= [A  X (0„  V 5),  4>n  A ((On  V B)  X (On  V 5))] 

= [A,B]. 


Remember  that  we  use  0n  and  On  to  denote  the  zero  and  unit  elements  in  the  belt  (Tn,  X,  V), 
respectively.  Let  [0„,G]  E P„.  Since  0n  x (G)+  = 0„  x (G')+  = 0„,  (0„,G)  ~ (0n,G')  for 
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any  G and  G' . So,  [0„,G]  is  unique  in  Pn.  We  denote  it  by  [0n,0„].  Also, 

[A,P]v[0n,0n]  = [A  X {0n  V 0n)  V (0n  V P)  X 0n,  <t>n  A (On  V 5)  X (On  V 0)] 

= 5] 

= [0„,  0n]  V [A,  B]. 

Thus,  [0n,0n]  is  the  identity  for  the  operation  V. 

Next,  we  need  to  show  that  (Pn,  x)  is  a semi-group  with  identity.  Note  that 

([A,B]x[C,D])x[E,F] 

= [A  x C,<t>n  K ((On  vB)x  (On  V D))\  x [P,P] 

= [A  X C X E,  (j)n  A (On  V ( (f)n  A (On  V B)  X (On  V T>))  X (On  V P))] 

= [A  x C x P,  <)>n  A (On  V P)  x (On  vD)x  (On  v F)] 

= [A,B}x([C,D}x(E,F]). 

Thus,  the  operation  X is  associative.  Note  also  that  [O„,0n]  is  in  Pn  and 

[A,  B]  X [On,  0„]  = [A  X On,  <t>n  A (On  V B)  X (On  v 0n)] 

= [A,B] 

= [On,0n]  x [A,B]. 

Therefore,  [O„,0n]  is  the  identity  for  the  multiplication  x. 

Finally,  we  show  that  the  operation  X is  distributive  over  the  operation  V. 

(A,B]x([C,D]s/[E,F\) 

= [A,  B]  X [C  X (On  v F)  V (On  VD)xUA  ((On  v D)  X (On  V F))] 

= [Ax  (C  X (On  v F)  V (On  VP)  X £),<£„  A ((On  VP)  X (On  V D)  X (On  V P))]. 

Note  that 


A X (On  V B)  X [cfrn  A ((On  V C)  X (On  V P)  X (0n  V P))]  + 
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= 4x(0„V5)x[4A(0„VCvBv(Cx5)V(flx5)v(CxBx  5))]+ 

= A X (On  V B)  X [C  V B V (C  X B)  V (B  X B)  V (C  X B X B)]+ 

= Ax(OnVB)x[CvBv(Cx  B)]+ 

= A x [C  V B V (C  x B)]+ 

= A x [<t>n  A ((On  V C)  x (On  V 5))]+. 

Thus, 

[A  X (On  V B),(f>n  A {(On  V C)  X (On  V B)  X (On  V 5))] 

= [A,  4>n  A ((On  V C)  X (On  V 5))]. 

Therefore, 

[A,5]x  [C,D]v[A,B]x  [E,F] 

= [AxC,<]>  /\  (On  V B)x  (On  V D)\  V [A  X E,(f)A  (On  V B)  X (On  V F)] 

= [A  X C X (On  v B)  X (On  Vf)V  (0n  V 5)  X (0n  V D)  X A X E, 

A (On  V B)x  (On  v D)  x (0n  V 5)  X (On  V F)] 

= [Ax  C x (On  vf)v  (On  vD)xAxE,4>/\  (0n  V5)  x (0n  v F>)  x (0n  V F)] 

= [A,B]x([C,D]v[E,F]). 

Similarly,  it  is  easy  to  verify  the  right  distributive  law. 

Q.E.D. 

Now,  let  T^x  denote  the  set  of  causal  recursive  templates  in  (R^,  rY00)y.  Define  a 
mapping  $ : T^-x  ~ * Pn  by  $(t)  = (\P(t^),  $(tx)).  Also,  define  a relation  ~ on  Txx  as 
follows: 


s ~ t 


if  and  only  if  $(s)  *(t) 
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Since  ~ is  an  equivalence  relation,  the  relation  ~ is  an  equivalence  relation  on  Txx-  We 
denote  the  set  of  all  equivalence  classes  over  T^x  by  T^x- 

Define  a mapping  <l>  : T^x  Pn  by  4>([t])  = [^(t^ ),  ^(t^)].  The  following  theorem 
characterizes  the  algebraic  structure  of  a class  of  non-linear  causal  recursive  templates. 

Theorem  4.23  4>  is  a belt  isomorphism  o/( T^Xi  0^,  V^)  onto  ( Pn , X,  V). 

Proof: 

If  [s],[t]  £ TxX,  then 
*([«]V-<[t])  = $([sv^t]) 

= 4([(s/I2I(0  V t^)  V (0  V A (0  V sx)0(O  V t^))]) 

= [®(s^0(O  V t^)  V (0  V s^)l3t^),  9(<(>  A (0  V 8^)121(0  V t«<))]. 

By  Theorem  4.20,  we  have 

k[*]v«  w) 

= [*(8/)  X (On  V *(t„))  V (On  V *(8^))  X 

<f>n  A (On  V $(8^))  X (On  V $(t^))] 

= [»(8^),  «(SX)]  V ®(t^e)] 

= *([8])V$([t]). 

Similarly,  we  have 

$([s]0^[t])  = $([s0^t]) 

= $([(s/[Y]t*,<M(0Vs^)[V](0Vtx))]) 

= [*(s/E]t*),¥(^A(0Vsx)E](0Vt.<))] 

= [*(8/)  X A (On  V ¥(sy))  X (On  V $(t^))] 

= [^(s^),  9(s^>]  X [®(t^),«(t^)] 

= ^([8])  x $([t]). 
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If  [A,  B]  6 Pn,  then  B is  strictly  upper  triangular  block  Toeplitz  with  upper  triangular 
Toeplitz  blocks.  We  know  that  $ is  1-1  and  onto  by  Theorem  4.20.  It  follows  that  t = 
($-1(A),  defines  a causal  recursive  template  in  Txx>  where  4>(t)  = [A,  B],  and  so 

4>  is  onto. 

Suppose  that  <l([s])  = 4>([t]).  Then,  4>(s)  ~ $(4).  By  definition,  we  have  s ~ t,  which 
implies  that  [s]  = [t].  Thus,  $ is  1-1  and  onto. 

The  unit  recursive  template  in  T^x*  denoted  by  0X,  is  given  by  = (0,0),  where  0 is 
the  empty  template  in  R*^.  By  definition,  $([0X])  = [O„,0„],  which  is  the  unit  element 
in  Pn. 


Q.E.D. 

Next,  we  show  some  basic  properties  of  the  recursive  operation  0.  From  the  above 
theorem,  we  immediately  have 


r0^(s0^t)  = (r0xs)0xt 


s0^t  = t0^s 

r0^(s  V^j  t)  = (r0xs)  (r0^t) 

for  any  r,s,t£  T^x- 

The  following  theorem  is  given  in  Zimmermann  [41]. 

Theorem  4.24  Let  X be  an  1 x n vector  of  indeterminates,  A an  1 X n vector  with  values 
in  R_oo,  and  B £ Mn.  If  B^  exists,  then  the  following  recursive  equation 

X = AM  (X  x B) 


has  the  unique  solution 


X = Ax  B+ 
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Theorem  4.25  If  a£  and  t 6 (R?oo,R?oo  )x  a causal  recursive  template,  then 

^(a[V]-(t)  = i/(a)  X 'fr(t7<)  X 

Proof: 

Let  b = alVJ-'t.  We  have 

Hb)  = (i/(a)  x ^(t/))  V (i/(b)  x 

Since  t is  causal,  is  upper  triangular  and  has  — oo  along  its  main  diagonal,  and  so 

^(t^)  is  stable.  Thus,  by  theorem  4.24,  we  have 

"(b)  = x $(t/)  x *(tx)+. 

Q.E.D. 

Lemma  4.26  If  A,  B E.  Tn  such  that  the  main  diagonal  elements  of  A and  B are  — oo,  then 

( <t>n  A ((0„  V A)  X (On  V B)))+  = A+xB+. 

Proof: 

If  A,  5 6 TB,  then  (A  V B)n  = V"=0  ^ x £n~\  and  thus  (A  V B)+  = Vo<,+j<„-i  A*  X BK 
Note  that  A‘  x B J = 0n  if  j + j > n.  It  follows  that 

A+  x B+  = \J  A*  x Bj 

°<»+i<n-i 

= (AV5)+. 

Therefore, 

(</>„  A ((On  V A)  X (On  V B)))+  = (<j>n  A (On  V A V B V (A  X B)))+ 

= (AvBv(Ax  B))+ 

= (A  V B)+ 


A+  x B+. 
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Q.E.D. 

Theorem  4.27  If  a G and  s,  t G (R?co,R^  )x  causal  recursive  templates,  then 

a0^(s0^t)  = (a[V]-<s)[V]^t. 

Proof: 

By  definition,  s|V]^t  = (s^0t/,(^  A (0  V s^)0(OVt^)).  By  Theorem  4.25,  Theorem  4.20, 

and  Lemma  4.26 

KaElx(s[2M))  = v(&)  X ®(s^|V]t^)  X [W(<£a  (0  V s^)0(O  V t^))]+ 

= v(&)  X 9(a^)  x ¥(t*)  X [<f>n  A ((On  V \P(s^))  X (On  V $(t^)))]+ 

= v( a)  X *(s,<)  X $(t/)  X (*(s^))+  X («(t^))+ 

= i/(aEUs)  x ¥(t*)  x (\P(t^))+ 

= i/((a|V]-<s)|Vl-<t) 

Q.E.D. 

Theorem  4.28  If  a,  b G R^  and  t G (R?^,  RX00)X  causal  recursive  template,  then 

(a  V b)EM)  = (aEM)  V (bBxt). 

Proof: 

By  Theorem  4.25, 

K(a  V b)|V|-<t)  = v(aVb)  X $(Q)  X f(t^)+ 

= (Ka)vHb))x%)xf(Q)+ 

= i/(a[V]^t  V bB_<t). 
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Q.E.D. 

Interestingly,  if  a 6 R?^,  and  s,t  G (R^,,  then  the  following  distributive  law 

a|V]-;(sV_<  t)  = (a0xs)  V (aH^t) 

does  not  hold  in  general.  For  example,  let  X be  an  1 X 3 point  set.  Let  a G RX  be  an  image 
and  s,t  G (R?^, RX00)X  two  simplified  recursive  templates  given  below 


a = 

1 

1 

i 

s = 

2 

t = 

3 

w 

Let  r = s[V]^t.  We  have 


r ^ 

T~ 

= 

4 

rx  = 

5 

3 

It  is  easy  to  compute  that 


a[V)^(s  V t) 


(a(V]^s)  V (a|V]-<t)  = 


58 


Thus,  a[V]-<(s  t)  ^ (a0^s)  V (a0xt). 

As  we  have  mentioned  in  Section  4.1,  all  the  results  here  would  be  valid  for  any  class  of 
quarter- plane  templates  as  well  as  causal  templates. 

Note  that  throughout  this  section,  we  assume  that  the  belt  operation  X is  -f.  Thus,  the 
results  are  about  the  recursive  operation  |V]X.  The  same  results  would  hold  for  the  recursive 
operation  if  we  change  the  belt  operation  x to  *. 

4.3  Max-Min  Transforms  and  Polynomials 

Discrete  linear  systems  are  usually  described  by  difference  equations,  which  then  may 
be  transformed  by  z-transform  into  polynomials  for  further  analysis.  In  this  section,  we  will 
apply  this  idea  to  morphological  filters.  We  show  that  man-min  transforms  can  be  defined, 
in  parallel  to  z-transform  in  linear  systems,  to  transform  morphological  filters  into  max-min 
polynomials.  Then,  in  Chapter  5,  we  will  develop  template  decomposition  techniques  based 
on  the  factorization  of  max-min  polynomials. 

Definition  Let  X C Z be  an  1 x n point  set,  and  a £ R^,.  The  max-transform  of  a, 
denoted  by  A4{a(i)}  = A(x),  is  defined  by 

A(x)  = V (a(0  ~ ix )• 

t=0 

The  min-transform  of  a G R^,  denoted  by  A4*{a(i)}  = A*(x),  is  defined  by 

A*(x)  = f\ (a(t')  - ix). 

t=0 

The  max-min  transform  has  many  similar  properties  that  the  usual  z-transform  has.  We 
will  prove  here  the  linearity  and  shifting  properties. 

Linearity: 

If  ki,  k2  £ R,  a,  b £ R^,  and  c,  d £ R^,  then 


M{(k!  + a(*))  V (*2  + b(i))}  = (fci  + A4{a(0})V(Jfc2  + A4{b(j)}) 
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Ai*{(ki  + c(i))  A {k2  + d(i))}  = (fca  + At*{c(i)})  A (jfc2  + A4*{d(i)}). 


Proof: 

Note  that  (R_oc,  V,-f)  is  a division  belt.  So, 

n— 1 

M{(k\  + a(i))  V (&2  + b(i'))}  = \J  [(&!  + a(i))  V (k2  + b(z))]  - ix 

i—0 

71—1 

= V [(*i + a(0  — ix)  V (k 2 + b(i)  — ix)] 

t=0 

7i—l  n— 1 

= V (*i  + a(0  " ix)  v V (^2  + b(i)  - ix) 

*'= 0 i=0 

= [*h  + V (a(0  ” ix)\  v [fe2  + V (mo  “ ix)l 

t'=0  j=0 

= (&i  + A(z))  V (k2  + P(a:)). 


Q.E.D. 


Shifting  Property: 

If  A:  > 0,  a G R-^oo,  and  b G R^,  then 

■Ad{a(i  — &)}  < At{a(i)}  — kx 

M*{b(i-k)}>M*{b(i)}-kx. 

Proof: 

71—1 

A(x)  — kx  = ( (a(i)  — ix))  — fcx 

t'=0 

71—1 

= \J  (a (i)  - ( i + k)x). 
i= o 

Substituting  j for  i + k and  noticing  that  a(i)  = — oo  for  i < 0,  we  have 

n+k—l 

A(x)  - kx  = \f  (a (j  - k)  - jx) 
j=k 
n+k— 1 

= V (aC?"  - k)  ~ ix) 

3=  0 
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n— 1 

> V (aM  ~k)~  ix ) 

i=o 

= A4{a(j  - fc)}. 

Q.E.D. 

The  proofs  for  the  linearity  and  shifting  property  of  the  min-transform  M*  are  similar. 
Difference  equations  are  used  to  describe  discrete  linear  systems  and  to  realize  recursive 
digital  filters.  We  will  use  max-min  equations  to  describe  recursive  morphological  filters. 
Definition  Let  X be  an  1 x m point  set,  and  a,  b £ R**,.  An  nth-order  dilation  filter  is 
defined  by 

MO  = V (ak  + a0  “ £))  V \J  ( bk  + b(i  - k)) 

k= 0 fc= 1 

and  if  a,  b £ R^,,  then  an  nth-order  erosion  filter  is  defined  by 

n— 1 n— 1 

MO  = A (a*  + a0  - M) A A (bk  + b(*  - *))> 

k-0  k=  1 

where  a*,  and  bk  are  filter  coefficients  in  R. 

Definition  The  system  function  of  a dilation  (or  erosion)  filter  is  defined  by 

H (x)  = B(x ) — A(x), 

where  A(z)  and  B(x)  are  the  max  (or  min)  transform  of  the  input  and  output  images, 
respectively. 

Lemma  4.29  The  inequality  x < a V x T b has  the  following  solution 

x < a - (0  V 6), 

where  a,b  £ R_oo>  and  x is  the  indeterminate. 

Proof: 

Let  x'  = a — (0  V b)  T r.  Substitute  x'  for  x in  both  sides  of  the  given  inequality.  Then, 


61 


1.  for  b > 0, 


Ihs  = a — b + r<a-\-r 
rhs  = flV(flfr). 

So,  Ihs  < rhs  is  always  true,  regardless  of  the  values  of  a and  r. 

2.  for  b < 0,  the  left  and  right-hand  side  of  the  given  inequality  will  be 

Ihs  = a + r 

rhs  = a V (a  + r + b)  = a + (0  V (r  + b)). 

Now,  if  r > 0,  then  (0  V (r  + 6))  < r,  and  whence 

rhs  = a + (0  V (r  + 6))  < a + r = Ihs. 

That  is,  the  inequality  does  not  hold.  On  the  other  hand,  if  r < 0,  then  clearly, 
Ihs  < rhs  for  any  a. 

Thus,  for  any  r < 0,  x = a — (0  V b)  -f  r is  a solution  of  the  given  inequality.  Therefore, 
^ < a — (0  V b)  is  the  solution  for  the  given  inequality. 


Q.E.D. 

Theorem  4.30  Given  an  nth-order  dilation  filter, 

n— 1 n— 1 

M0  = V + a(* - M) v V (M  + M*  - &)) 

k= 0 fc=l 

its  system  function  is  of  the  form 

n— 1 n— 1 

H{x)  = B(x)  - A{x)  < \J  (ak  - kx ) - (0  V \/  (bk  - kx)). 

k= 0 k= 1 
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Proof: 

Applying  the  max-transform  and  using  the  linearity  and  shifting  property,  we  get 

B{x)  = M{\f(ak  + a(i-  k))  V \/  (bk  + b(t  - k))} 
k= 0 k= 1 

n— 1 n— 1 

= V (ak  + Ad{a(i  - A:)})  V \J  (bk  + M{b(i  - k )}) 

k= o jt=l 

n— 1 n— 1 

< V ( ak  + A(»)  - kx)  V \/  + _0(x)  - kx) 

k= 0 fc=l 

n— 1 n— 1 

= (A(x)  + \J  (ak-  kx))  V ( B(x ) + \f  (bk  - kx)). 

k= 0 k=l 

Solve  this  inequality  for  B(x),  using  Lemma  4.29,  we  have 

B(x)  < A(x)  + \J  (ak  - kx)  - (0  V \/  (bk  - kx)). 

k= 0 fc=l 

Q.E.D. 

Note  that  the  system  function  is  the  difference  of  two  max-polynomials,  which  charac- 
terizes the  given  recursive  morphological  filter. 

Similar  conclusions  can  be  made  about  min-transform. 

Lemma  4.31  The  inequality  x > a A (a:  + b)  has  the  following  solution 

x > a — (0  A 6), 

where  a,b  E R+co,  and  x is  the  indeterminate. 

Theorem  4.32  Given  an  nth-order  erosion  filter, 

n—1  n-1 

MO  = A (°*  + M*  ~ M) A A (bk  + M*  “ M) 

fc=0  k= 1 

its  system  function  is  of  the  form 

H*(x)  = Y*(x)  - X*(x)  > A {ak  - kx)  - (0  A f \(bk - kx)). 

k= 0 k= 1 


CHAPTER  5 

TEMPLATE  DECOMPOSITION 


Both  linear  convolution  and  morphological  methods  are  widely  used  in  image  processing. 
One  of  the  common  characteristics  among  them  is  that  they  both  require  applying  a tem- 
plate to  a given  image,  pixel  by  pixel,  to  yield  a new  image.  In  the  case  of  convolution,  the 
template  is  usually  called  convolution  window  or  mask;  while  in  mathematical  morphology, 
it  is  referred  to  as  structuring  element.  Templates  can  vary  from  one  to  another  in  their 
weights,  sizes  and  shapes,  depending  on  the  specific  applications. 

Intuitively,  the  problem  of  template  decomposition  is  that  given  a template  t,  find  a 
sequence  of  smaller  templates  1 1 , . . . , tn  such  that  applying  t to  an  image  is  equivalent  to 
applying  tj, . . . , tn  sequentially  to  the  image.  In  other  words,  t can  be  algebraicly  expressed 
in  terms  of  ti, . . .,tn. 

One  of  the  concerns  with  template  decomposition  is  that  some  of  the  current  image 
processors  can  only  handle  very  small  templates  at  a time.  For  example,  ERIM’s  Cyto- 
Computer  [20]  cannot  deal  with  templates  of  size  larger  than  3 X 3 on  each  pipeline  stage. 
Thus,  a large  template  has  to  be  decomposed  into  a sequence  of  3 x 3 or  smaller  templates 
before  it  can  be  applied  to  an  image  on  CytoComputer. 

A more  important  motivation  for  template  decomposition  is  to  speed  up  template  op- 
erations. For  large  convolution  mask,  the  computation  cost  resulting  from  implementation 
can  be  prohibitive.  However,  in  many  instances,  this  cost  can  be  significantly  reduced  by 
decomposing  the  masks  or  templates  into  a sequence  of  smaller  templates.  For  instance,  the 
linear  convolution  of  an  image  with  an  n x n template  requires  n2  multiplications  and  n2  — 1 
additions  to  compute  a new  image  pixel  value;  while  the  same  convolution  computed  with 
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<te 

an  1 X n row  template  followed  by  an  n X 1 column  template  takes  only  2 n multiplications 
and  2(n  — 1)  additions  for  each  new  image  pixel  value.  This  cost  saving  may  still  hold 
for  parallel  architectures  such  as  mesh. connected  array  processors  [42],  where  the  cost  of  a 
convolution  is  proportional  to  the  size  of  the  template. 

The  problem  of  template  decomposition  has  been  investigated  by  several  people.  Ritter 
and  Gader  [3]  presented  some  very  efficient  methods  for  decomposing  DFT  and  general 
linear  convolution  templates.  Wiejak  and  Buxton  [43]  proposed  a method  to  decompose 
a 2-D  (or  higher  dimension)  Marr-Hildreth  convolution  operator  into  two  1-D  convolution 
operators.  Zhuang  and  Haralick  [44]  gave  an  algorithm  based  on  tree  search  that  can 
find  an  optimal  2-point  set  decomposition  of  a morphological  structuring  element  if  such  a 
decomposition  exists. 

In  this  chapter,  we  will  first  discuss  the  issues  of  template  decomposition  for  both  re- 
cursive and  nonrecursive  templates.  We  then  give  a necessary  and  sufficient  condition  for 
decomposing  separable  rectangular  templates,  and  a simple  algorithm  to  find  such  decom- 
positions. Next,  we  will  develop  methods  for  decomposing  symmetric  convex  templates  and 
spherical  templates,  and  show  that  the  computation  cost  of  the  resulting  decomposition  can 
be  reduced  by  an  order  of  magnitude. 

5.1  Fundamentals  of  Template  Decompositions 

The  following  properties  of  the  template  operations  in  the  image  algebra  [3,45]  provide 
a mathematical  basis  for  template  decomposition. 

If  a G Rx,  and  r,  s,  t G (RX)X,  then 

a©(s©t)  = (a©s)©t 
a©(s  + t)  = (affis)  + (a©t) 


r©(s  + t)  = (r©s)  -f  (r©t). 
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If  a G RX.,  and  r,  s,t  6 (R?w)^)  then 

a0(s0t)  = (a0s)0t 

a0(s  V t)  = (a0s)  V (a0t) 
r0(s  V t)  = (r|V]s)  V (r0t). 

These  properties  are  important  because  they  help  to  explore  the  possibilities  of  comput- 
ing template  operations  in  different  ways.  As  far  as  the  efficiency  is  concerned,  the  goal  of 
template  decomposition  is  to  find  the  most  efficient  way  to  implement  a template  operation 
with  a given  template.  For  instance,  if  we  know  that  r = s©t,  then,  by  the  associative  law, 
we  could  apply  s and  t sequencially  to  a instead  of  computing  a©r  directly  since  s and  t 
are  much  smaller  than  r in  general. 

In  Chapter  4,  we  proved  some  similar  associative  and  distributive  laws  involving  re- 
cursive templates  and  operations.  They  can  be  used  for  the  decomposition  of  recursive 
templates. 

Suppose  a G RX,  and  r,  s,  t G (Rx,  Rx)x.  We  have  proved  the  following  in  Chapter  4, 

a®x(s©xt)  = (affixs)©xt 

a©x(s  +x  t)  = (a©xs)  + (a©xt) 
r©x(s  +x  t)  = (rffixs)  +x  (r©xt). 

If  a G R?oo,  and  r,s,t  G (R?^,  RX00)X,  then 

a[V]x(s[V)xt)  = (a0xs)|V]xt 
r0_<(s  Vx  t)  = (r0xs)  Vx  (r|V]xt). 

Definition  A ffi- decomposition  of  a template  t G (Rx)x  is  a sequence  of  templates  1 1 , . . ., 
tn  G (Rx)x,  such  that  t = ti  © . . . ©tn. 
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Definition  A 0- decomposition  of  a template  t G is  a sequence  of  templates 

ti,...,tn  G (R^00)x,  such  that  t = ti0  . . . 0tn. 

By  using  both  associative  and  distributive  laws,  a weak  version  of  template  decomposi- 
tion can  be  defined  also. 

Definition  A weak  © -decomposition  of  a template  t G (RX)X  is  a sequence  of  templates 
ti , • • • , tfcj , . . . , tfcn_1+i , . . . , tfcn  G (RX)X  such  that 

t = (ti  © • • • ©tfcl ) + (t*1+1  © . . . ©tfc2 ) + ...+  (tkn_1+1  © . . . ©tfcj. 

For  example,  if  a G RX,  and  t G (Rx)x  have  a weak  decomposition  t = ti©t2+t3©t4, 
then  we  can  compute  affit  as  follows 

a©t  = ((a©ti)ffit2)  + ((a©t3)©t4). 

In  general,  the  decomposition  of  a template  t,  if  possible,  is  prefered  to  its  weak  decom- 
position because  usually  more  time  and  space  would  be  involved  in  computing  and  applying 
a weak  decomposition. 

The  decomposition  of  recursive  templates  are  similarly  defined.  For  example,  a ©_-- 
decomposition  of  a recursive  template  t G (RX,RX)X  is  a sequence  of  templates  tx, ...  ,tn  G 
(Rx)x,  such  that  t = ti  ©x  . . . 

With  the  concept  of  template  decomposition  defined,  we  show  next  how  to  decompose 
some  commonly  used  templates. 

5.2  Decomposition  of  Separable  Templates 

Definition  An  invariant  template  t G (FX)Y  of  finite  support  is  called  a rectangular  tem- 
plate if  its  support  «S(ty)  (or  S-o^ty)  is  of  rectangular  shape  at  each  target  pixel  y G Y.  A 
recursive  invariant  template  t G (Fx,  FY)Y  of  finite  support  is  called  a recursive  rectangular 
template  if  both  and  1 — (or  (0  V t^))  are  rectangular  templates. 
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Rectangular  templates  are  the  simplest  and  most  commonly  used  templates  in  image 
processing.  Ideally,  to  speed  up  template  operations,  we  would  like  to  decompose  an  arbi- 
trary rectangular  template  into  two  1-dimensional  templates,  namely  a row  template  and  a 
column  template.  Thus,  for  a given  m X n rectangular  template,  the  number  of  arithmetic 
operations  required  for  each  template  operation  on  each  pixel  can  be  reduced  from  0(mn ) 
to  0(m  + n ). 

Definition  A rectangular  template  is  called  separable  if  it  can  be  expressed  as  the  com- 
position of  a row  template  and  a column  template. 

Let  r G (FW)Y  be  an  m x n rectangular  template  with  weights  W{j,  i = 1, . . .,m;  j = 
1,  • • • , n.  We  say  that  r is  0- separable  if  there  exist  an  1 X n row  template  s G (Fw)x  and  an 
TO  X 1 column  template  t G (FX)Y  such  that  r = s0t;  ^separable  if  r = s|V]t.  Similarly, 
if  r G (Fx,  Fx)x  is  a recursive  rectangular  template,  then  We  say  that  r is  0 ^-separable  if 
there  exist  recursive  templates  s,  t G (FX,FX)X,  where  and  are  row  templates  and 
t^  and  t,<  are  column  templates,  such  that  r = s0xt;  M^-separable  if  r = s0xt. 

Next,  we  give  a sufficient  and  necessary  condition  for  the  separability  of  rectangular 
templates. 

Theorem  5.1  If  r E (RXoc)Y  is  an  m X n rectangular  template  with  weights  rij,  where 
i = l,...,m;  j = 1, ...,»,  then  r is  ^-separable  if  and  only  if  for  all  1 < i < to,  and 
1 < j <n, 

rij  ~ rij  = r.i  - ru  (5.1) 

Proof: 

Suppose  that  r is  0-separable  and  r = s[V]t,  where  s is  an  1 X n row  template  with  weights 
si, . . .,5n  and  t is  an  to  X 1 column  template  with  weights  ti,...,tm.  By  the  definition  of 
template  composition,  r,y  = f;  + Sj,  for  i = 1, . . . , to,  and  j = 1, . . . , n.  Thus, 


rij  rlj  — ( ti  T sj)  (tl  Tsj) 
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= U-h 

— (U  + si)  — (h  + si) 

= rn  - rn. 

Now,  assume  that  Equation  (5.1)  holds.  Define  an  1 X n template  s and  an  m X 1 
template  t as  follows: 


sj  ~ ri  j (5.2) 

<«  = rn  - rn  (5.3) 

Thus,  for  j = 1, . . .,  n and  i — 1, . . .,  m 

rij  = (rn  - rn)  + rxj 

— ti  T Sj . 

Therefore,  r = s[V]t,  and  r is  [VJ-separable  as  desired. 

Q.E.D. 

This  yields  a straight  forward  method  for  testing  and  decomposing  0-separable  tem- 
plates. Given  anmxn  rectangular  template  r,  it  takes  (to-  l)n  additions  and  (TO-l)(n- 1) 
comparisons,  according  to  the  equation  (5.1),  to  see  whether  r is  separable  or  not.  If  it  is 
separable,  then  one  can  easily  construct  the  corresponding  row  and  column  templates  by 
the  equations  (5.2)  and  (5.3)  as  given  above. 

Theorem  5.2  If  r € (RX)Y  is  a to  X n rectangular  template  with  weights  r,y,  where  i = 
1, . . . , to;  j = l,...,n,  then  r is  © -separable  if  and  only  if  for  all  1 < i,  1 < to,  and 
1 < j,  1 < n, 


rij/rlj  = Tn/m 


(5.4) 
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The  proof  is  similar  to  that  of  Theorem  5.1. 

Notice  that  the  condition  (5.4)  is  equivalent  to  saying  that  the  rank  of  (r,y)  is  1. 

Sometimes,  especially  when  defining  large  templates,  it  is  rather  convenient  to  define  a 
rectangular  template  of  size  mxn  by  a function  W(x , y ) of  two  variables  over  an  m X n grid. 
We  call  such  function  a weight  function.  On  the  other  hand,  any  real  valued  function  of 
two  variables  defined  on  a finite  rectangular  grid  defines  a real  valued  rectangular  template, 
whose  support  is  the  domain  of  the  function.  Thus  the  separability  of  a rectangular  template 
can  be  reduced  to  the  separability  of  its  weight  function.  The  following  results  are  obvious. 


Theorem  5.3  If  t is  a real  valued  rectangular  template  of  size  mxn,  and  Wt(x,y)  its 
weight  function,  where  x = 1, . . . , m and  y = 1, . . . , n,  then  t is  © -separable  if  and  only  if 
Wt{x,y)  = f(x)  • g(y)  for  some  real  valued  functions  f and  g;  t is  M-separable  if  and  only 
if  Wt(x,  y)  = f(x)  + g(y)  for  some  real  valued  functions  f and  g. 

Example  We  define  a (2m  + 1)  X (2n  + 1)  paraboloid  template  t by  the  function 

W(x,y)  = k(x2  + 1 12) 


where  k is  a constant,  -m  < x < m,  and  -n  < y < n.  Then  t is  El-separable. 
Example  An  n X n Marr- Hildreth  template  t is  defined  by 


1 x2  + y'v_(x 
7T( 7*  ' 2<72 


1—fl  _ . e— y2/2cr2  , _i -x2/2<t2(,  _ f_)e-y>/2*2 

r <t4  ' „2>e  +27TCT4  1 


27 rcr4 


= Mx)f2(y)  + f3(x)f4(y) 


where  — n < x < n,  and  — n < y < n. 

Clearly  t is  not  separable  but  it  still  has  a very  efficient  weak  ©-decomposition,  which 
is  the  sum  of  two  separable  templates,  as  follows: 


t = (ti  ©t2)  + (t3  ©t4) 
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where  tj  is  defined  by  /,  for  i = 1,2, 3, 4. 

Finally,  since  a recursive  template  actually  consists  of  two  nonrecursive  templates,  its 
separability  can  be  derived  from  the  separability  of  the  two  nonrecursive  templates. 

Theorem  5.4  If  r E (FX,FX)X  is  a quarter-plane  rectangular  recursive  template,  then  r 
is  (V]  ^-separable  if  and  only  if  both  and  0 V are  [V]  - separable ; r is  © ^-separable  if 
and  only  if  both  r ^ and  1 — are  © -separable. 

Proof: 

We  will  prove  the  ©^-separable  case.  Suppose  that  r is  ©_< -separable.  Then,  there 
exist  recursive  templates  s,t  £ (Fx,  Fx)x,  where  and  are  row  templates  and  and 
are  column  templates,  such  that  r = s©^t.  By  definition, 

r/  = s*©^ 

= 1 - (1  - sx)©(l  - t^j) 

So, 

1 - ^ = (1  “ sx)©(l  - tx) 

Note  that  (1  - s<)  is  still  a row  template  and  (1  - tx)  a column  template.  Thus,  both 
and  1 — are  ©-separable. 

Now,  if  both  and  1 — are  ©-separable,  then  there  exist  row  templates  si,s2  and 
column  templates  ti,t2  such  that 


rj  = si©tx 


1 - rx  = s2©t2 


So, 


r^  = 1 - s2©t2 
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We  define  two  recursive  templates  s and  t by 

s = (si,l-s2) 
t = (t1?l-t2) 

Since  r is  a recursive  template  and  1 — = s2ffit2,  s and  t are  well  defined  recursive 

templates.  It  is  easy  to  see  that 

s©^t  = r 

Hence,  t is  ©^-separable. 

Q.E.D. 

5.3  Decomposition  of  Symmetric  Convex  Templates 

Symmetric  convex  templates  are  often  used  in  morphological  filters.  In  this  section,  we 
show  how  to  decompose  these  symmetric  convex  templates  with  respect  to  the  El  operation. 

The  concept  of  convex  set  is  usually  defined  in  n-dimensional  Euclidean  space  Rn.  In 
image  processing,  an  image  has  to  be  digitized  on  a subset  of  Zn.  Hence,  it  is  not  clear 
what  a digital  convex  set  should  be.  In  this  section,  we  first  define  what  a digital  convex 
set  is  in  Z2  and  what  a convex  template  is.  We  then  define  a set  of  special  symmetric 
convex  templates,  called  cross  templates,  and  show  that  all  boolean  symmetric  convex 
templates  can  be  decomposed  as  a sequence  of  boolean  cross  templates.  We  can  actually 
show  a stronger  result,  that  is,  a boolean  symmetric  template  is  convex  if  and  only  if  it 
is  the  composition  of  a sequence  of  boolean  cross  templates.  Here  again  we  only  consider 
translation  invariant  templates. 

Let  X be  a subset  of  Z2.  We  say  that  a point  x = (a;i,x2)  G X is  4- adjacent  to  a point 
y = (j/i>  2/2)  G X if  | xx  — y1  | -f  | X2  — 3/2  |=  1.  A point  x G X is  4-connected  to  a point 
y £ X if  there  is  a path  x = xo,xi, . . .,xn  = y such  that  xj  6 X and  xj  is  4-adjacent  to 
xi-l  for  1 < i < n.  X is  4-connected  if  for  any  pair  of  points  x,y  G X,  x is  4-connected 
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to  y.  X is  symmetric  if  it  is  symmetric  with  respect  to  both  x-axis  and  y-axis.  X is  finite 
if  it  is  contained  in  a disk  of  finite  radius.  If  X is  finite  and  4-connected,  then  a point  x is 
called  an  extreme  point  of  X if  x is  in  X but  x does  not  lie  between  any  two  points  of  X. 
In  other  words,  x has  at  most  one  horizontal  neighbor  and  at  most  one  vertical  neighbor  in 
X. 

Definition  A subset  X of  Z2  is  called  convex  if  it  is  finite,  4-connected,  and  satisfies  the 
following  two  conditions: 

1.  The  set  of  extreme  points  of  X forms  a convex  polygon  H in  R2 

2.  For  any  x £ Z2,  x 6 X if  and  only  if  x G H 

Definition  An  invariant  template  t £ is  a convex  template  if  its  support  ^^(ty) 

is  convex  for  each  y£Y.  t is  called  a symmetric  convex  template  if  its  support  is  both 
symmetric  and  convex. 

It  is  clear  that  a convex  set  in  Z2  is  uniquely  defined  by  its  extreme  points,  thus  so  is 
the  support  of  a convex  template.  From  the  definition  of  convex  template,  we  immediately 
have  the  following 

Lemma  5.5  If  t is  a symmetric  convex  template  and  (xXt  yx ),...,  (xn,  yn)  the  extreme 
points  of  S — oq (ty ) in  the  first  quadrant,  where  xx  < . . . < xn  and  y = (0,0),  then  for 

1 = 2,...,n 

min{z,-  - y,-_!  - yi}  = 1 

and  for  i = 2, . . . , n — 1 

Vi- 1 ~ Vi  > yi  ~ Vi+ 1 

%i— 1 *£«+l 

Definition  A cross  template  cross (i,j)  is  defined  as  the  union  of  a row  template  of  size 

2 * i + 1 and  a column  template  of  size  2 * j + 1 on  their  target  pixels,  where  min{i,y}  is 


either  0 or  1. 
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Some  examples  of  cross  templates  are  shown  below: 


lii 


0 

0 

z'  N 
xV 

0 

0 

cross(l,0)  cross(l,l) 


0 

1 

1 

1 

1 

0 

cross(2, 1) 

Clearly  all  cross  templates  are  convex  templates.  Furthermore,  they  are  irreducible  in 
the  sense  that  they  cannot  be  decomposed  into  a sequence  of  smaller  templates  (except  the 
simple  rows  and  columns).  A far  more  important  property  is  that  all  boolean  symmetric 
convex  templates  can  be  represented  as  the  0-composition  of  a sequence  of  boolean  cross 
templates.  A template  t G (R?00)Y  is  called  boolean  if  all  of  the  weights  in  its  support  are 
0’s. 

Theorem  5.6  Any  boolean  symmetric  convex  template  t G (R?0Q)Y  has  a decomposition 
t = crossi  0 . . . 0crossm,  where  cross i(*i , Ji.), . . . , cross m(im,jm)  are  boolean  cross  tem- 
plates. 

Proof: 

We  prove  this  by  induction  on  n,  the  number  of  extreme  points  of  S-o^ty)  in  the  first 
quadrant.  If  n=l,  let  {x\,y{)  be  the  extreme  point  in  the  first  quadrant.  Then  t is  a 
boolean  rectangular  template  which  is  separable  and  t = cross^x,  0)0 cross(0,  y\).  For 
n > 2,  let  (xi,  yi), . . .,(xn,  yn ) be  the  extreme  points  of  5_oo(ty)  in  the  first  quadrant  such 
that  xi  < ...  < xn.  We  define  a boolean  cross  template  as  crossi(x2  — xi,f/i  — 2/2)-  It 
can  be  shown  by  using  Lemma  5.5  that  (t  El  crossi  )0  crossi  = t and  t E3  crossi  is  still  a 
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boolean  symmetric  convex  template,  whose  support  has  n - 1 extreme  points  (xi,y2),  (2:3  - 
Ax,  2/3), . . .,  (x„  — Ax,  ?/„),  where  Ax  = X2  — Xj.  Thus  by  induction  hypothesis,  tEjcrossj 
has  a 0-decomposition  of  a sequence  of  cross  templates.  Hence,  so  does  t. 

Q.E.D. 

If  t £ (RXqo)Y  }s  a boolean  symmetric  convex  template  and  cross(i,y)  is  an  arbitrary 
boolean  cross  template,  then  t0cross (i,j)  is  still  a boolean  symmetric  convex  template.  So, 
the  0-composition  of  any  finite  number  of  boolean  cross  templates  is  a boolean  symmetric 
convex  template. 

Theorem  5.7  A template  is  a boolean  symmetric  convex  template  if  and  only  if  it  can  be 
represented  as  the  0 -composition  of  a sequence  of  cross  templates. 

The  set  of  all  boolean  symmetric  convex  templates  from  Z2  to  Z2  are  generated  by  the 
set  of  boolean  cross  templates  and  closed  under  the  operations  0 and  [A]. 

The  following  algorithm  is  derived  from  the  Theorem  5.6  and  can  be  used  to  decompose 
an  arbitrary  boolean  symmetric  convex  template  into  a sequence  of  cross  templates. 

Algorithm  1 (Decomposition  of  Boolean  Symmetric  Convex  Templates) 

Input:  A boolean  symmetric  convex  template  t £ (R^00)Y. 

Output:  A sequence  of  cross  templates  cross1(i1,y1), . . . , crossn(i„,  jn)  such  that 

t = crossi(ii,y1)0  . . . 0cross„(in,  jn) 

Step  1:  Generating  the  boundary  of  the  template 

Let  r be  the  radius  of  the  given  template  along  the  x-axis.  We  need  only  to  store  the  first 
quadrant  of  the  boundary  due  to  the  symmetry.  Let  c[0],  c[l], . . . , c[r]  be  an  array  of  integers, 
representing  the  boundary  in  the  first  quadrant. 
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Step  2:  Computing  the  supports  of  the  cross  templates 

From  the  array  c,  build  an  n X 2 index  array  idx  such  that  each  entry  in  idx  defines  the 
support  of  a cross  template. 

n = 0 

pre_i  = 0 

for  i = 0 to  r 

if  (c[i]  > c[i+l])  then 
n = n + 1 

idx[n,l]  = c[pre_i]  - c[i] 
idx[n,2]  = i - pre_i 
pre_i  = i 
end  if 
end  for 
n = n + 1 
idx[n,l]  = c[r] 
idx[n,2]  = 0 

Step  3:  Generating  the  sequence  of  cross  templates 

Once  the  index  array  idx  is  set  up  as  above,  the  cross  templates  we  are  looking  for  will  be 
crossi(icfr[l,  1],  ida;[l,2]), . . cross n(idx[n,  1],  idx[n,  2]) 


where  the  weights  of  each  cross  template  are  zeros. 
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Example  The  followings  are  an  elliptical  disk  template  with  major  axis  21  and  minor  axis 
15,  and  a disk  template  of  radius  21.  They  are  generated  by  composing  a sequence  of  cross 
templates  obtained  by  the  above  algorithm. 


Figure  1.  Examples  of  symmetric  convex  templates 


5.4  Decomposition  of  Spherical  Templates 

In  many  morphological  filters,  a disk  template  or  spherical  template  is  desirable  because 
it  defines  an  equal-distance  neighborhood. 

Definition  An  invariant  template  t E (R?00)Y  is  called  a disk  template  if  it  is  boolean 
and  its  support  <S_oo(ty)  is  a digital  disk. 

Definition  An  invariant  template  t E (R?00)Y  is  called  a spherical  template  if  its  support 
(S_oo(ty)  is  a digital  disk  and  the  weights  in  its  support  define  a digital  half  sphere. 

Now  the  problem  is  that  given  an  integer  r,  we  need  to  find  a sequence  of  templates 
1 1 , . . . , tk  such  that  the  composition  tj  0 . . . 0tjc  gives  rise  to  a disk  (or  spherical)  template 
of  radius  r.  Note  that  a disk  template  is  a special  boolean  symmetric  convex  template,  and 
thus  can  be  decomposed  by  the  algorithm  given  in  the  last  section. 
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We  show  next  that  a spherical  template  can  also  be  decomposed  as  a sequence  of  cross 
templates  with  various  weights.  The  algorithm  for  decomposing  spherical  templates  is 
quite  similar  to  the  previous  algorithm  except  that  we  have  to  consider  how  to  compute  the 
weights  of  each  cross  template  generated. 

Algorithm  2 (Decomposition  of  Spherical  Templates) 

Input:  A spherical  template  t 6 (RXOQ)Y  of  radius  r. 

Output:  A sequence  of  cross  templates  crossi(ij,  ji), . . . , crossn(in,  jn)  such  that 

t = cross! (*i,  Ji)EI . . . Elcrossri(tn,;n) 

Step  1:  Generating  a digital  circle 

Use  Bresenham’s  circle  algorithm  [46]  to  generate  an  optimal  digital  circle  of  radius  r.  We 
need  only  to  store  the  first  quadrant  of  the  circle  due  to  the  symmetry.  Let  c[0],  c[l], . . . , c[r] 
be  the  array  of  numbers,  representing  the  circle  in  the  first  quadrant.  Note  that  c[0]  is  equal 
to  r. 

Step  2:  Computing  the  supports  of  the  cross  templates 

From  the  array  c,  build  an  n X 2 index  array  idx  such  that  each  entry  in  idx  defines  the 
support  of  a cross  template. 

n = 0 

pre_i  = 0 

for  i = 0 to  r 

if  (c[i]  > c[i+l])  then 
n = n + 1 

idx[n,l]  = c[pre_i]  - c[i] 
idx[n,2]  = i - pre_i 
pre_i  = i 
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end  if 
end  for 
n = n + 1 
idx[n,l]  = c[r] 
idx[n,2]  = 0 

Step  3:  Generating  the  sequence  of  cross  templates 

Once  the  index  array  idx  is  set  up  as  above,  the  cross  templates  we  are  looking  for  will  be 

crossi(zda;[l,  1],  idx[  1,2]), . . . , crossn(idx[n,  1 \,idx[n,  2]) 

where  the  weights  of  each  cross  template  crossk(idx[k,  l],  idx[k,  2])  can  be  computed  as  fol- 
lows: 

for  i = -idx[k,l]  to  idx[k,l] 

weight_cross_k[i,0]  = idx[k,i]  - idx[k,2] 
end  for 

for  j = -idx [k,2]  to  idx[k,2] 
weight_cross_k(0, j)  = 0 
end  for 

weight_cross_k(0,0)  = idx[k,l] 
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Example  The  following  is  a spherical  template  of  radius  16,  which  is  generated  by  com- 
posing a sequence  of  cross  templates  obtained  by  the  above  algorithm. 


Figure  2.  An  example  of  spherical  template 

Next,  we  show  that  the  decompositions  of  disk  and  spherical  templates  can  result  in 
very  efficient  template  operations. 

When  a template  is  applied  to  an  image,  the  amount  of  computations  involved  is  pro- 
portional to  the  size  of  the  support  of  the  template.  Thus,  the  larger  the  template  is,  the 
more  computations  would  be  required. 

Theorem  5.8  Let  t = crossi(ti,  j\)M  . . . Bicross  n(in,jn)  be  a decomposition  of  some  disk 
(or  spherical)  template  t 6 (R^0O)Y  of  radius  r.  IfT(r)  = card(S- oo(t))>  the  size  of  the 
support  oft,  and  To(r)  = card(<S_oo(crossi))  -(-...  + card(<S_oo(crossn)),  then  Tr>(r)  < 
5r  + 1. 

Proof: 

Note  that  n < r + 1 for  any  r > 1,  and 

*i  + *2  + • • • + in  — r 


jl  + J2  + • • • + 7n  = r 
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It  follows  that 

Td{t)  = card(5_00(crossi))  + . . . + card(5_00(crossn)) 

= (2(«i  + jl)  + 1)  + ■ • • + (2  (in  + jn)  + 1) 

= 2(ii  + *2  + . . . + in)  + 2(il  + J2  + • • • + jn)  + n 

= 4r  + n 
< 5r  + 1 

Q.E.D. 

It  is  easy  to  see  that  T(r ) > nr2.  Thus  after  the  decomposition,  the  template  operation 
can  be  faster  by  an  order  of  magnitude. 

5.5  Max-polynomials  and  Template  Decomposition 

In  this  section,  we  will  briefly  review  the  max-polynomials  and  show  how  they  can  be 
used  to  decompose  morphological  templates. 

Definition  A max-polynomial  P(x)  in  indeterminate  x is  defined  by 

n 

P(x)  = V (a«  + rix)  - (ao  + rox)  V (ai  + r\x)  V • • • V (an  + rnx ) 

«'=o 

where  ao, . . . , an  are  in  R_oo,  r0, . . . , rn  are  in  Z,  and  tq  < • • ■ < rn. 

We  will  denote  the  set  of  all  max-polynomials  in  x with  coefficients  in  R_oo  by  R^^fa:]. 
Definition  If  P(x)  = a0V(aa  + a;)V-  • -V(am  + mx)  and  Q(x)  = 60V(6i  + x)  V- • -\/(bn  + nx) 
are  both  in  and  m < n,  then 

P(x)  V Q(x)  = c0  V (ci  + a:)  V • • • V (cn  + nx ) 

where  c;  = ai  V 6,-,  for  0 < i < n.  We  assume  that  a,-  = — oo  for  m < i < n. 

Definition  If  P(x)  = a0V (a i + x)V  - • 'V(am+mi)  and  Q(x)  = 60V(6!  + *)V-  • -V(bn  + nx) 
are  both  in  R^jx],  and  m < n,  then 

P(x)  + Q(x)  = c0  V (ci  + x)  V • • • V (cm+n  -I-  (rn  + n)x) 
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where  c;  = (a,-  + b0)  V (a,_i  + h)  V • • • V (a0  + &,),  for  0 < i < m + n. 

Definition  Let  P(x)  = a0  V + *)  V • • • V (an  + nx)  be  in  R*^®].  If  P(x)  ± -oo  and 
an  ^ — oo,  then  the  degree  of  P(x)  is  n. 

Let  P(x)  = V"=o(ai  + rix)  be  a max-polynomial.  We  define  fa,  and  fa  as  follows: 
fa(x)  = ak  + rkx  (k  = 0, . . . , n) 

If  n > 1: 

n 

fa{x)  = \J  fa(x)  (k  = 0,  ...,n) 

i= 0,t/fc 

ak- 1 — ak  ,,  x 

Pk  = («  = 1, . . n) 

rk  ~ rk- 1 

Thus,  P(x)  can  be  written  as 

n 

F(x)  = V ^*(*) 

i=0 

One  significant  difference  between  a usual  polynomial  and  a max-polynomial  is  that  a 
max-polynomial  may  have  redundant  terms,  which  may  be  dropped  without  affecting  the 
functionality  of  the  max-polynomial.  For  example,  since  l + x<2V(l  + 2x),Vx  £ R, 

2 V (1  + x)  V (1  + 2x)  = 2 V (1  + 2x) 

Vx  £ R. 

Definition  Let  P(x)  = Vr=o(ai  + rix ) be  a max-polynomial.  If  fa  and  fa  are  defined  as 
above,  then  fa  is  inessential  in  P{x)  for  a particular  i (0  < i < n)  if  the  following  holds  for 
all  x: 

<f>i(x)  < Mx) 

Otherwise  we  will  say  that  ^ is  essential. 

Since  r0  < r\  < ■ ■ • < rn,  when  x — > -oo  and  y —*■  -foo, 


fa(x)  ~ fa{x)  = (a0  - a,)  + (r0  - r,)x  ->  +oo 


for  i = 1, . . .,n 
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4>n(x)  - = (an  - a,)  + (rn  - rt)x  -*■  +oo  for  i = 0, . . n - 1. 

Thus,  4> o and  (f>n  are  always  essential. 

A max-polynomial  P(x)  is  an  irredundant  expression  if  n = 0 or  if  n > 1 and  <f>k  is 
essential  in  P(x)  for  all  0 < k < n. 

Note  that  if  </>,(x)  is  inessential,  then  d>,(x)  may  be  removed  from  the  definition  of  P(x) 
without  affecting  the  function  value  of  P(x),  i.e.,  we  may  write  P(x)  = ipi(x).  Here,  we 
have  to  distinguish  P(x)  as  a formal  expression  and  as  a function  of  x.  Thus,  by  removing 
<l>i(x),  P(x)  would  be  changed  as  a formal  expression  although  its  functionality  may  be  the 
same. 

The  following  theorem  is  from  Guninghame-Green  and  Meijer  [47]. 

Theorem  5.9  A max-polynomial  P(x ) is  an  irredundant  expression  if  and  only  if  either 
n < 1,  or  n > 1 and  /?i  <•••</?„. 

5.6  Factorization  of  Max-polynomials 

A max-polynomial  may  be  factored  into  the  sum  of  linear  terms.  For  notational  conve- 
nience, we  will  let  Factor(P(x ))  and  Expand(P(x))  denote  the  factorization  of  P(x)  and 
the  expansion  of  P(x),  respectively.  For  instance,  if  P(x)  = (1  V x)  -f  (2  V x),  then 

Expand(P(x ))  = 3 V (1  + x)  V (2  + x)  V 2x 
= 3 V (2  + x)  V 2x. 

The  following  theorem  proved  by  Cuninghame-Green  and  Meijer  [47]  is  an  analogue  of 
the  fundamental  theorem  of  algebra,  in  which  /3,’s  play  the  roles  of  roots  of  a polynomial. 

Theorem  5.10  Let  P(x)  he  an  irredundant  max-polynomial. 

n 

p{x)  = Y(a,-  + r;x) 

t=0 
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Then  P(x ) has  a unique  linear  factorization 

n 

Factor(P(x ))  = an  + r0x  + ]^(r,  - r;_i)(ft  V x) 

i= 1 
n 

= an  + r0x  + ^((at_a  - a,)  V (r,-  - 
1=1 

Note  that  if  we  view  a max-polynomial  P(x)  as  a function,  then  we  may  drop  any 
inessential  terms  to  make  it  irredundant.  So,  functionally,  every  max-polynomial  P(x)  has 
a unique  linear  factorization. 

For  the  purpose  of  template  decomposition,  we  are  interested  in  max-polynomials  as 
formal  expressions  instead  of  their  functionalities.  We  know  that  morphological  template 
with  respect  to  [V)  can  be  represented  as  a max-polynomial,  and  thus  may  be  decomposed 
by  the  factorization  of  the  corresponding  max-polynomial.  The  indeterminate  x here  simply 
plays  the  role  of  a place  holder.  The  factorization  corresponds  to  the  template  decomposition 
whereas  the  expansion  corresponds  to  the  template  composition. 

However,  due  to  the  irredundancy  of  max-polynomials,  the  above  theorem  may  not  be 
used  directly  for  template  decomposition,  which  can  be  illustrated  by  the  following  counter- 
example. 

Example  Let  P(x)  = 4 V (3  + x)  V (1  + 3z).  Then, 

ft  = 4-3  = 1,  ft  = (3  - l)/2  = 1. 

So,  ft  < ft,  and  hence  P(x)  is  irredundant.  By  the  factorization  theorem,  we  have 

Factor(P(x ))  = 1 + (1  Vi)  + (2  V 2x). 

But  if  we  expand  the  resulting  factorized  expression,  we  will  get 

Expand(F actor (P(x)))  = 4 V (3  + x)  V (2  + 2z)  V (1  + 3x), 

which  is  not  identical  to  the  original  P(x)  as  a formal  expression  although  they  have  the 
same  functional  values. 
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The  problem  here  is  that  the  expansion  of  the  factorized  expression  of  P(x ) may  intro- 
duce some  extra  inessential  terms  that  may  not  exist  in  the  given  max-polynomial  P(x ) in 
the  first  place.  Remember  that  the  indeterminate  x is  simply  a place  holder  for  the  weights 
of  a template.  So,  if  we  think  P(x)  as  a template  and  its  factorization  as  the  decomposition 
of  the  template,  then  this  decomposition  would  be  incorrect  because  we  cannot  get  the  given 
template  back  by  convolving  the  decomposed  templates.  Fortunately,  we  may  get  around 
this  problem  by  adding  some  more  restrictions  on  P(x). 

Theorem  5.11  Let  P(x)  be  an  irredundant  max-polynomial 

n 

p{x ) = V (a,-  + r,x) 

»'= o 

such  that  ai  ^ — oo  for  all  i = 0, . . . , n,  and  r4-  — r8_i  = c for  all  i = 1, . . . , n,  where  c is 
some  constant  integer.  Then  P(x)  has  a unique  linear  factorization 

n 

Factor(P(x ))  = an  + r0x  + ^(r,-  - r,_i)(/?j  V x) 

t'=i 

n 

= an  + r0x  + ]T((at_i  - a,)  V (r;  - r,_i)x) 

«=i 

Furthermore, 

n 

Expand(F  actor  (P(x)))  = P(x)  = ( a,-  + r8x) 

i=o 

Proof: 

The  first  part  is  a special  case  of  the  factorization  theorem.  We  prove  the  second  part  by 

induction  on  n.  It  is  a trivial  case  for  n = 0.  Assume  that  P(x ) = Vi=o(a«  + rix ) an^ 

k 

Expand(F actor (P(x)))  = P(x)  = -f  r,x) 

»=o 

for  all  k < n.  Then,  let  P(x)  = V”=o (ai  + r{x).  Expand  its  factorization  and  use  the 
induction  hypothesis,  we  have 


Expand(Factor(P(x ))) 
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n 

— Expand(an  + tqx  + ^((a;_i  — a,)  V cx )) 

1 = 1 

n— 1 

= Expand(an-\  + r0x  + [^((a.-i  - a;)  V cx)]  + (an  - a„_i)  + ((a„_i  - an)  V cx)) 

i— 1 

n— 1 

= [ V (a*  + r«*)]  + (°  V («n  - «n-l ) + cx) 

i=0 

n— 1 n— 1 

= [ V (a«  + r«a:)] v [ V (°*  + (a«  - a«-i)  + (r«  + c)x)] 

i=0  t=0 

n— 1 n— 2 

= [ V (°*  + r«*)]  V [ V (a*  + (fln  - an-l)  + (ri  + c)x)]  V («n  + rnX). 

1=0  i= 0 

Since  < @2  ■■•<  /3n  and  rt-  — = c for  i = 1, . . . , n,  we  have 


ao  ^1  ^ ^1  ^ — 1 


Hence, 

®t— 1 T (®n  ®n— 1 ) <'- 

for  all  i = 1, . . . , n — 1.  It  follows  that 

n— 1 n— 2 n— 1 

[ V (a*  + rix)\ v [ V (a«  + (a«  ~ a«-i)  + (r*  + c)*)]  = [ V (a*  + r«,a;)]- 

j=l  i= 0 t=l 

Thus, 

n 

Expand(Factor(P(x)))  = \J  (a,  + r,a;)  = P(x). 

i= l 

Q.E.D. 

The  above  theorem  gives  a sufficient  condition  for  decomposition  of  1-dimensional  mor- 
phological templates  but  it  is  not  a necessary  condition.  For  example,  let  P(x)  = 4 V (3  + 
2x)v(2-|-3x)v(l-i-5x).  Clearly,  P(x)  is  redundant,  and  thus  we  cannot  apply  the  above 
theorem.  However,  it  is  easy  to  check  that 

Expand^  1 T (1  V 2x)  + (2  V 3a;))  = P(x). 

Here,  (2  + 3a:)  is  an  inessential  term  in  P(x). 

The  general  procedure  to  decompose  an  1-dimensional  template  t is  as  follows: 
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1.  Find  the  max-polynomial  representation  P(x ) of  t. 

2.  If  P(x)  satisfies  the  conditions  of  Theorem  5.11,  then  get  the  factorization  and  convert 
all  terms  back  to  templates,  which  will  be  the  decomposition  for  t;  otherwise,  goto 
next  step. 

3.  Remove  all  inessential  terms  in  P(x ) and  get  the  factorization  of  the  irredundant 
expression  P'(x ) of  P(x). 

4.  If  Expand(Factor(P\x )))  = P(x),  then  convert  all  terms  back  to  templates,  which 
will  be  the  decomposition  for  t;  otherwise,  t cannot  be  decomposed  into  a set  of 
two-point  templates. 

Note  that  we  may  combine  the  method  presented  here  with  the  techniques  given  in 
Section  5.3.  For  example,  given  a paraboloid  template  r below, 


0 

3 
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3 

0 

3 
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7 

6 

3 

4 

7 
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7 
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3 

6 

7 

6 

3 

0 

3 

4 

3 

0 

We  know,  from  Section  5.3,  that  it  is  separable.  In  particular,  we  have  r = s[V]t,  where 
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s = 


liJ 


t = 


o 


3 

3 


0 


Next,  let  P\(x)  and  P^x)  be  the  max-polynomial  representations  of  s and  t,  respectively. 

Pi(x)  = (0-2x)v(3-x)V4V(3  + x)v(0-|-2x) 

P2(x)  = Pi(x). 

By  Theorem  5.11, 


Factor(P\(x ))  = -2x  + (-3  V x)  + (-1  V x)  + (1  V x)  + (3  V x) 

= (—3  — x V 0)  + ( — 1 — x V 0)  + (1  V x)  + (3  V x). 


Thus, 


r > 

-3 

0 

El 

-i 

o 

El 

i j 

0 

El 

3 

0 

t = 

-3 

El 

-l 

El 

a 

0 

0 
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From  Section  4.3,  we  know  that  a recursive  morphological  template  may  be  represented 
as  the  difference  of  two  max- polynomials,  thus,  may  be  decomposed  by  the  method  pre- 
sented here.  Also,  the  factorization  techniques  for  max- polynomials  can  be  applied  to 
min-polynomials  as  well. 

We  conclude  this  chapter  by  pointing  out  that  template  decomposition  is  not  only 
necessary  in  the  case  where  special  image  processing  hardware  cannot  handle  large  templates 
at  a time,  but  also  desirable  when  the  problem  of  efficient  computation  is  concerned  in  real- 
time applications.  In  this  chapter,  we  have  shown  that  general  template  decomposition  is 
possible  based  on  the  associative  properties  of  template  operations  in  the  image  algebra. 
If  both  associative  and  distributive  properties  are  used,  a weak  decomposition  may  be 
derived  also.  The  necessary  and  sufficient  conditions  as  well  as  algorithms  are  given  for 
decomposing  rectangular  shaped  templates  and  symmetric  convex  templates.  We  have 
presented  efficient  methods  for  decomposing  disk  and  spherical  templates  which  are  often 
used  in  mathematical  morphology.  The  amount  of  computation  time  required  to  apply  such 
templates  can  be  reduced  by  at  least  an  order  of  magnitude.  Finally,  we  have  shown  that 
template  decomposition  may  be  obtained  by  factorization  of  max-polynomials. 


CHAPTER  6 

APPLICATIONS  IN  SIGNAL  AND  IMAGE  PROCESSING 


In  this  chapter,  we  present  some  applications  of  the  recursive  template  operations  in 
signal  and  image  processing.  We  show  how  easily  the  image  algebra  can  be  used  to  express 
various  recursive  algorithms  in  image  processing. 

6.1  Moving  average  filter 

A fc-order  causal  moving  average  filter  is  defined  by: 

^ k— 1 

n)=  -rJ2x(n-i 

K t=0 

where  x(n ) and  y(n ) are  the  input  and  output  sequences,  respectively. 

This  filter  is  a non-recursive  filter.  But  we  may  turn  it  into  an  equivalent  recursive  filter 
as  follows: 

y(n)  = y(n  - 1)  + ^(x(n)  - x(n  - k)) 

Note  that  each  output  needs  k additions  for  the  non-recursive  version  of  the  average 
filter,  whereas  only  2 additions  per  output  are  required  for  the  recursive  version  no  matter 
how  large  the  k is. 

Let  X be  an  1 X n point  set.  Let  a,b  £ RX,  be  the  input  and  output  images,  respectively. 
Define  a recursive  template  t G (RX,  Rx)x  as  follows: 


— 1/k 

k 

u = 

i 
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Then,  the  image  algebra  expression  for  the  moving  average  filter  would  be 

b :=  a©^t 

Since  two-dimensional  average  filter  is  separable,  it  can  be  realized  by  the  cascade  of  two 
one-dimensional  average  filters,  one  along  the  vertical  direction  and  one  along  the  horizontal 
direction.  For  a,  k X k causal  average  filter  with  k > 2,  a recursive  implementation  would 
be  more  efficient  on  a sequential  machine. 

6.2  Recursive  Filters 

Various  linear  recursive  filters,  such  as  low-pass,  high-pass,  etc.,  are  used  in  image 
processing  [31].  A 2-D  recursive  filter  is  characterized  by  a 2-D  difference  equation 

y(m,n)=  ^ a(k,l)x(m  — k,n  — l)  — ^ b(k,l)y(m  — k,n  — l) 

k,leRa  k,l£Rb-(  0,0) 

where  x(m,  n ) and  j/(m,  n ) are  input  and  output  sequences,  Ra  and  Rb  are  the  regions  of 
support  for  constants  a(k,  l ) and  b(k,  l ),  respectively. 

If  we  define  a recursive  template  t in  such  a way  that  the  weights  of  and  t_/  are  given 
by  the  filter  constants  a(k,l ) and  b(k,  /),  respectively,  then  a general  recursive  filter  can  be 
specified  in  the  image  algebra  by 

c©^t 

where  c is  assumed  to  be  the  input  image. 

So  the  template  operations  © and  ©^  are  sufficient  to  describe  all  linear  recursive 
and  nonrecursive  filters.  Furthermore,  they  are  more  general  than  recursive  filters  since 
templates  can  be  variant  in  the  image  algebra  and  nonlinear  template  operations  may  be 
used  to  express  nonlinear  filters. 
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6.3  Predictive  Compression 

When  a large  image  a is  transmitted  over  a communication  link  of  limited  capacity,  it  is 
desirable  to  compress  the  image  first.  One  easy  and  fast  way  to  do  this  is  to  use  predictive 
compression  technique  (DPCM)  [32,40]. 

Let  X C Z2  and  a G R.  Suppose  that  a is  correlated.  Then  it  is  possible  to  estimate 
or  predict  a pixel  value  a (i,j)  in  terms  of  its  previous  neighbors.  For  example, 

a (*,  j)  = c\&(i  - 1,  j)  + c2a(i,  j - 1) 

where  a (i,j)  is  the  predicted  value  for  a (i,j). 

Define  the  differential  image  e G RX 

e(*,i)  = a(i,j)-a(t,j) 


Note  that  e contains  less  information  (hence  requires  less  storage)  because  images  generally 
contain  areas  of  almost  constant  or  slowly  changed  gray  levels,  and  in  these  regions,  pixel 
values  of  e are  close  to  zero. 

Now  define  a template  s G (RX)X  and  a simplified  recursive  template  t G (RX,RX)X 
below 


-.5 

-.5 

^ SN 

.5 

.5 

Then  the  compression  process  can  be  described  by  the  expression 

e = /(a©s) 


where  / is  some  nonlinear  quantizer. 
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The  original  image  a can  be  reconstructed  approximately  by  the  receiver  from  the  quan- 
tized differential  image  e as  follows 

a — e©^t 

where  -<  is  the  forward  scanning  (from  left  to  right  and  top  to  bottom)  order  defined  on 
X.  Note  that  actually,  t = s-1.  Thus,  the  information  may  only  be  lost  in  the  quantization 
process. 

6.4  Distance  and  Medial  Axis  Transforms 

The  distance  transform  is  a nonlinear  transform,  that  may  be  implemented  either  re- 
cursively or  nonrecursively  [48].  They  were  first  proposed  by  Rosenfeld  and  Pfalz  [36,49]. 
The  recursive  version  works  much  faster  than  the  nonrecursive  version  if  only  sequential 
machines  are  used. 

Let  a be  a binary  image  with  +oo  for  feature  pixels  and  0 for  non-feature  pixels.  The 
distance  transform  of  a is  a gray  level  image  b such  that  each  pixel  value  b (t,j)  is  the 
distance  between  a (i,j)  and  the  nearest  non-feature  pixels.  The  distance  can  be  measured 
as  Euclidean  distance,  or  city  block,  chess  block,  etc.,  subject  to  the  specific  applications. 

The  distance  transform  can  be  used  to  compute  the  medial  axis  transform  [36].  Let 
X C Z2,a  G {0,+oo}x,  and  b E RX. 

First,  the  distance  transform  b of  a is  given  by 

b = (a[A]  ) O tr 

where  -<  is  the  forward  scanning  (from  left  to  right  and  top  to  bottom)  order  defined  on  X 
and  >-  is  the  backward  scanning  (from  right  to  left  and  bottom  to  top)  order  defined  on  X. 

The  weighted  medial  axis  W(a)  of  a is  the  set  of  local  maxima  of  b,  which  can  be 
computed  by 


W(a)  = b * x=fl(b0s  — b) 
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Finally,  the  reconstruction  of  b from  W(a)  can  be  done  by 

b = (W(a)|V|^(-t))[V]^t* 

Here,  the  templates  s,t  may  be  defined  according  to  the  specific  distance  measure  being 
used  [48].  For  example,  if  we  want  to  measure  the  city- block  distance,  then  we  may  define 
the  template  s and  the  simplified  recursive  template  t as  follows: 


1 

1 

6.5  Smoothing  and  Edge  Detection  Using  Recursive  Filtering 

Multiscale  smoothing  and  detecting  edges  of  an  image  can  be  done  efficiently  by  using 
the  recursive  filtering  technique  [50], 

Let  X be  an  1 X n point  set  and  a 6 Rx  the  input  image.  Then,  the  smoothed  image 
b e Rx  is  given  by 

b(i)  = b1(i)  + b2(i)  for  i = l,..., n 

where 

bj(*)  = k(a.(i ) + e~a(a  — l)a(j  — 1))  + 2e~ah\(i  — 1)  — e~2ohi{i  — 2) 
b2(i)  = k(e~a(a  + l)a(i  + 1)  - e~2“a (i  + 2))  + 2e-ab2(*  + 1)  - e~2ab 2(z  + 2) 

and 

J:_  (l~e-«)2 

1 -(-  2ae-“  - e~2a 

Note  that  a is  a scaling  factor.  By  adjusting  the  value  of  a,  one  can  effectively  control  the 
size  of  the  smoothing  operator  and  thus  the  amount  of  noise  suppression  without  increasing 
the  amount  of  computation  for  each  output  pixel. 


94 


Now,  if  we  define  two  recursive  templates  s and  t as  follows 


7^s 

s*  = 

•S2 

= 

s4 

S3 

\ y 

7^s 

h 

II 

Y 

■*> 

f > 

U 

where 

si  = k,si  = ke~a(a  — 1),S3  = 2e_a,S4  = —e~2a 

and 

ti  = ke~a(a  + 1),  = -e~2o,,t3  — 2e~a  ,t4  — —e~2a 

Then,  the  image  algebra  expression  for  the  smoothing  operator  will  be 

bi  :=  a©xs 
b2  :=  a©>_t 
b :=  bi  + b2 

Smoothing  of  2-D  images  can  be  achieved  by  cascade  of  two  1-D  smoothing  operators, 
one  along  the  horizontal  direction  and  one  along  the  vertical  direction.  Specifically,  let 
X is  an  n X n point  set.  Let  the  recursive  template  s,t  as  before,  and  let  sv,tv  be  two 
corresponding  vertical  recursive  templates.  Then,  the  smoothing  of  a 2-D  image  a £ R is 
as  follows 


b 


i 


a©^s 
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b2  :=  a©^t 
b3  :=  bi  + b2 
bi  :=  b3ffi^sv 
b2  :=  b3ffi^tv 
b :=  bi  + b2 

The  same  recursive  filter  structure  can  be  used  to  realize  a multiscale  first  derivative 
operator  If  we  change  the  weights  of  the  recursive  templates  s and  t as  follows 

5i  = 0,52  = 1,  s3  = 2e-cv,  s4  = — e~2a 

and 

= l,/2  = 0,  t3  = 2e~a , <4  = —e~2a 
6,6  Sequential  Labeling 

One  fundamental  task  in  image  processing  is  that  given  a binary  image  of  four-connected 
blobs,  assign  each  blob  a different  label  so  that  further  processing  can  be  made  on  each 
individual  blob  separately.  Among  many  other  methods,  the  sequential  labeling  algorithm 
[33]  is  found  to  be  the  fastest  one  on  conventional  computers. 

The  idea  of  this  algorithm  is  to  scan  the  given  image  recursively  from  left  to  right  and  top 
to  bottom  with  an  L-shaped  template.  Let  X £ Z2  be  an  m x n point  set.  Let  a £ {0, 1}X 
and  b £ Rx  be  the  input  image  and  output  image,  respectively.  Then,  the  image  algebra 
formulation  of  the  sequential  labeling  algorithm  is  given  by 

b :=  a©^t 

where  the  t is  a variant  parameterized  recursive  template  defined  as  follows: 


if  a (2/1?  j/2)  = 1 then 
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if  a(t/i  - 1, y2)  = 1 and  a (yi,y2  - 1)  = 0 then 

Mtfi  - 1,2/2)  = 1 

elseif  a(j/i  — 1,  y2)  = 0 and  a(r/i,  y2  - 1)  = 1 then 
Mj/i>2/2  - 1)  = 1 

elseif  a(j/i  — 1,  y2)  = 1 and  a(t/1?  y2  — 1)  = 1 then 
- l,lfe)  = 1 

b(t/!  - 1,  y2)  is  equivalent  to  y2  - 1) 

else 

cur. color  = cur. color  + 1 
, 2/2)  = cur. color 

endif 

endif 

where  cur. color  is  assumed  to  be  a global  variable  initialized  to  0. 

A simple  post  processing  with  lookup  table  can  be  applied  to  change  all  the  equivalent 
labels  to  a unique  label. 

6.7  Dithering  with  Error  Diffusion 

Digital  halftoning  or  spatial  dithering  is  the  method  of  rendering  the  illusion  of  con- 
tinuous tone  pictures  on  display  devices  that  are  capable  of  producing  only  binary  picture 
elements.  Dithering  with  error  diffusion  is  a popular  dithering  algorithm  [34].  It  is  a recur- 
sive algorithm.  The  input  image  is  scanned  in  forward  raster  scanning  order.  A threshold 
is  fixed  at  0.5  assuming  that  the  gray  values  in  the  input  image  vary  from  0 to  1.  The 
output  pixel  value  is  either  0 or  1,  depending  whether  the  input  pixel  value  is  greater  than 
the  threshold  or  not.  The  difference  between  the  thresholded  output  pixel  value  and  the 
input  pixel  value  is  referred  to  as  the  error,  which  will  then  be  propagated  to  a weighted 
neighborhood  along  the  scanning  direction.  These  past  errors  will  be  added  to  the  future 
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input  pixel  values  to  give  a correction  factor  so  that  errors  may  be  diffused  and  a smooth 
dithering  output  image  may  be  obtained. 

The  error  diffusion  algorithm  can  be  formulated  in  the  image  algebra  by  generalized 
recursive  template  operation  as  follows.  Let  X C Z2  be  an  m X n point  set,  a G Rx  the 
input  image,  and  b G {0, 1}X  the  output  image.  We  define  a simplified  recursive  template 
t G (RX,RX)X  as  follows 


1/16 

5/16 

3/16 

7/16 

iM 

Let  7 = + and  define  two  binary  functions  oj  and  o2  by 


°i(z,2/)  = xy 


and 

°2 {x,y) 

Then  the  image  algebra  expression 


( x — 1 )y  if  x > 0.5 
xy  otherwise 


b :=  x>o.5(a®xt) 

gives  a formal  specification  of  the  error  diffusion  dithering  algorithm. 

Note  that  if  we  let  c = a®xt,  then  for  each  y G X, 

c(y)  = H a(x)o1t7<y(x)+  Y1  c0)°2W(z) 
xe5(t^y)  zes(txy) 

= J2  a(x)t7<y(x)  + 5Z  (c(z)  - X>o.5(c(z))txy(z) 
xe-s(t^y)  ze-s(t^y) 

The  image  c can  be  viewed  as  the  intermediate  error-corrected  image,  which  will  then 


be  thresholded  to  generate  the  output  image. 


CHAPTER  7 

CONCLUSION  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


We  have  established  a theoretical  framework  based  on  the  image  algebra  for  studying 
recursive  operations  in  image  processing.  By  introducing  recursive  templates  and  recursive 
template  operations  into  the  current  body  of  the  image  algebra,  we  have  enhanced  the  ex- 
pressive power  of  the  image  algebra  and  shown  that  many  interesting  results  and  techniques 
can  be  derived  and  used  in  image  processing. 

In  particular,  we  have  shown  that 

1.  By  introducing  the  notion  of  recursive  templates  and  recursive  template  operations, 
the  image  algebra  can  be  used  to  express  many  useful  recursive  image  processing  al- 
gorithms in  a formal  and  yet  straight  forward  manner.  The  recursive  template  opera- 
tions are  capable  of  representing  both  linear  and  nonlinear  recursive  filters  commonly 
encountered  in  image  processing. 

2.  The  image  algebra  provides  an  algebraic  model  for  computations  of  recursive  opera- 
tions in  image  processing.  We  have  characterized  the  quarter-plane  linear  recursive 
filters  and  the  quarter-plane  nonlinear  recursive  filters  and  shown  that  they  are  iso- 
morphic to  the  algebraic  objects  ring  and  belt,  respectively. 

3.  The  computation  of  these  recursive  template  operations  are  generally  not  globally 
parallel  in  the  sense  that  the  computations  of  some  output  image  pixel  values  may 
depend  on  some  other  output  pixels.  But  in  many  cases,  partially  parallel  computation 
can  be  achieved. 
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4.  Research  in  developing  algorithm  for  factoring  max-polynomials  has  applications  to 
template  decomposition  in  parallel  image  processing.  Some  sufficient  conditions  are 
found  for  decomposing  recursive  morphological  templates  based  on  the  factorization 
of  max-polynomials. 

5.  Symmetric  convex  templates  are  commonly  used  in  mathematical  morphology,  which 
include  disk  templates  and  spherical  templates.  We  have  shown  that  the  class  of 
boolean  symmetric  convex  templates  can  be  generated  from  the  class  of  boolean  cross 
templates  and  are  closed  under  the  dilation  and  erosion  operations.  We  have  developed 
the  algorithms  to  decompose  an  arbitrary  boolean  symmetric  convex  template  into  a 
set  of  cross  templates  as  well  as  the  algorithm  to  decompose  spherical  templates. 

The  author  would  like  to  make  the  following  suggestions  for  further  research  and  would 
like  to  be  informed  of  any  progress  made  about  the  following  problems. 

The  first  question  is  to  generalize  the  results  in  this  dissertation.  In  Chapter  4,  We  have 
established  the  algebraic  structures  for  linear  quarter-plane  recursive  filters  and  nonlinear 
quarter-plane  recursive  filters.  We  know  from  Gader’s  dissertation  and  Davidson’s  disser- 
tation that  algebraic  structures  exist  for  any  linear  and  nonlinear  nonrecursive  filters.  Can 
algebraic  structures  be  found  for  non-symmetric  half-plane  recursive  filters  and  spatially 
variant  recursive  filters?  Since  these  filters  are  more  general  than  quarter-plane  recursive 
filters,  results  related  to  them  would  be  more  interesting  and  useful  in  image  processing. 

We  have  not  touched  the  existence  problem  of  local  decomposition  for  recursive  tem- 
plates in  this  research.  This  problem  has  been  solved  by  Gader  and  Davidson  for  linear 
and  nonlinear  nonrecursive  templates.  In  the  context  of  recursive  templates,  the  problem 
is  that  given  a recursive  template  and  a network  of  processors  interconnected  by  some 
communication  links,  decompose  the  template  into  a sequence  of  templates  such  that  each 
resulting  template  is  local  with  respect  to  the  processor  network.  Under  what  conditions 
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could  such  a local  decomposition  exist?  In  particular,  what  partial  orders  and  types  of 
processor  networks  would  permit  such  local  decomposition? 

Another  area  of  research  is  to  develop  more  general  template  decomposition  techniques. 
In  Chapter  5,  we  have  presented  techniques  for  decomposing  arbitrary  boolean  symmetric 
convex  templates.  Could  this  be  carried  out  to  all  gray-level  symmetric  convex  templates? 
We  have  only  solved  the  problem  for  spherical  templates,  which  are  a class  of  very  special 
gray-level  symmetric  convex  templates.  In  Chapter  5,  we  have  also  developed  sufficient 
conditions  and  techniques  for  decomposing  one-dimensional  recursive  templates  based  on 
the  factorization  of  max-polynomials.  In  the  light  of  the  fact  that  there  is  no  fundamental 
theorem  of  algebra  for  polynomials  of  two  or  more  variables,  it  is  unlikely  that  one  can  find 

a similar  theorem  for  max-polynomials  of  two  or  more  variables.  But  special  methods  may 

» 

be  developed  to  factor  some  special  max-polynomials  of  two  variables,  and  then  be  used  to 
decompose  two-dimensional  morphological  templates. 
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